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FOREWORD 

The SPACE-A program described in this report is a modified version 
of the trajectory and observation generation portion of the Sequential Position 
and Covariance Estimation (SPACE) program originally acquired from NASA. 
This document contains a detailed analysis of the equations and models used 
by the program, a brief description of the routines used in the program, as 
Well as instructions for its operation. The principal applications of SPACE-A 
are the prediction of space-vehicle trajectories and the generation of observa- 
tional data plus other trajectory related information. 

There have been numerous changes made in the original program in 
addition to its adaptation to the IBM 7030 and its subsequent debugging. The 
modified SPACE-A program is the result of the efforts of the authors and 
L. E. Wilkie. S. Schwartz aided in the debugging and rewriting of certain 
routines. Special thanks and appreciation are due to R. K. Squires, D. S. 
Woolston, and other members of the Special Projects Branch, Theoretical 
Division of the Goddard Space Flight Center from whom the SPACE program 
was obtained. 

The work reported in this document was performed by The MITRE 
Corporation, Bedford, Massachusetts, for the Directorate of Planning and 
Technology, Electronic Systems Division, of the Air Force Systems Command 
under Contract AF 19(628)-5165. 



REVIEW AND APPROVAL 



Publication of this technical report does not constitute Air Force approval of 
the report's findings or conclusions. It is published only for the exchange and 
stimulation of ideas. 

ANTHONY P. TRUNFIO 

Technical Advisor, Development Engineering Div. 

Directorate of Planning and Technology 
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ABSTRACT 

This report describes the SPACE-A program presently available for 
operational use, SPACE-A is the trajectory and observation generation 
portion of the Sequential Position and Covariance Estimation Program (SPACE) 
which is currently under development. The document contains a detailed 
analysis of the equations and models used by the program, a brief description 
of routines used in the program, as well as instructions for its operations. 
The principal applications of SPACE-A are the prediction of the space- 
vehicle trajectories and the generation of observational data plus other 
trajectory related information. 
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SECTION I 
OBJECTIVES 

INTRODUCTION 

The SPACE trajectory determination program was originally devel- 
oped for the Special Projects Branchy Theoretical Division, Goddard 
Space Flight Center by the Sperry Rand Systems Group. It was designed 
as a comprehensive trajectory determination and tracking program for 
both orbital and deep space flights. The original program consisted 
of three major modes of operation: 

(1) SPACE-A, designed for trajectory and observation generation 
without statistical processing. 

(2) SPACE-Bl, designed for statistical estimation of the six 
state variables describing the position and velocity at 
prescribed points of the trajectory. 

(3) SPACE-B2, designed for statistical estimation of the six 
state variables describing position and velocity plus up 
to 20 additional states selected from a group of vari- 
ables which permit estimation of dynamic biases (para- 
meters affecting orbital motion) or observational biases 
(affecting measurements). 

The version of the SPACE program received by MITRE contained 
numerous errors in programming. SPACE-Bl and SPACE-B2 had never been 
debugged and many of the routines of SPACE-A were unworkable. In ad- 
dition, the programming had to be made compatible with the IBM 7030. 
Therefore, a large effort was put into checking the theory, debugging 
the program, and revising or rewriting certain routines. Many check- 
out runs were necessary during and after the debugging of the program. 

Since the SPACE-B2 mode includes the capabilities of the SPACE-Bl 
mode, the SPACE-Bl mode was not debugged. The original SPACE-B2 
(which was not operational when received by MITRE) employed two dif- 
ferent statistical estimation techniques: 

(1) a minimum variance sequential filtering technique, the so- 
called "Kalman filter", and 



(2) a non-recursive batch processing technique. 

Because a number of trajectory estimation programs at MITRE 
already employ the batch processing technique, the second option of the 
SPACE-B2 mode was left in a non-operational status. The first option 
of SPACE-B2» which incorporates the Kalman flJter in the trajectory 
estimation procedure, is presently being debugged and checked out. 
The theory, programming, and results of its operation wil] be pub- 
lished in a succeeding document. The theoretical development and 
much of the programming of the trajectory generation portion of 
SPACE-B2 is identical with that used in SPACE-A. 

This document deals with the modified SPACE-A trajectory generci- 
tion program written in Fortran IV and is presently operational on 
the IBM 7030 at MITRE. This report is based on the original docuncnts'- ''■' ' J 
published by Sperry Rand for NASA. 

Section I contains a brief discussion of the objectives and capa- 
bilities of SPACE-A. Section II discusses the theoretical background, 
as well as the equations and methods employed by the program Section 
III consists of a user's guide for operating the program. Section fV 
gives the program structure and a brief description of the function 
of each subroutine. 

The larger part of the report is contained in Section II since 
understanding or following the programming is really a matter of fol- 
lowing the appropriate equations. Although parts of Section II arr 
original, the contents of pertinent sections of the SPACE Analytic 
Manual^ -^ have been freely used or modified to fit the description 
of the present program. Sections III and IV closely follow the format 
of the original User's Manual L^-l and Programmer's Manual *- -* ; however, 
since a number of changes in programming have been made, only Section 
III and Section IV of this report should be used in operating the modi- 
fied SPACE-A trajectory generation program. 



CAPABILITIES 

The mcxiified SPACE-A trajectory generation may be run in any of 
the following modes: 

(1) Trajectory generation - the computation of position and 
velocity of the vehicle at regular time points along a 
trajectory. 

(2) Observation generation - the computation of certain 
ground-based observations of the vehicle at regular 
time points. 

(3) Simulated data - the computation and writing onto 
tape of the observations in a format suitable for 
processing data by SPACE-B2. 

(4) Visibility computation - the time of initial view 
of the vehicle by a given station, the total time 
in sight by the station, and the observations 
while the vehicle is in sight. 

The required input consists of the initial date and initial time 
of a particular run and the initial position and velocity of the vehi- 
cle. When the effect of drag is to be included which is usually the 
case, the effective area and mass of the vehicle and the magnitude of 
solar flux must be specified. When the effect of direct solar radia- 
tion is included, an effective area pertaining to solar radiation must 
be specified. For observation calculations, the specification of the 
station location on the surface of the Earth is necessary. The de- 
tails of the units and format of these required inputs along with 
other optional input features are discussed in Section III. The types 
of output data, which consist primarily of trajectory and/or observa- 
tion information, is also given in Section III. 

The philosophy behind the modified SPACE-A trajectory generation 
program is that its primary use would be for Earth orbital trajec- 
tories. Toward this objective, its dynamic model includes the effects 
of the primary terra and oblateness perturbations of the Earth's grav- 
ity, the solar, lunar, and planetary gravitational perturbations, the 
effect of direct solar radiation and drag, as well as such dynamic ef- 
fects as precession, nutation, and the daily rotation of the Earth. 



The model used for drag is a combination of the U.S. Standard 
Atmosphere 1962 and the Harris-Priester model for the upper atmosphere, 
which depends on the magnitude of solar flux. The model for Earth 
gravitational oblateness may include up to nine zonal harmonics, four 
sectorial harmonics » and twelve additional tesseral harmonics. The 
types of ground station observations that may be specified are given 
in Section II (see OBSERVATIONS). 

Although the modified SPACE-A program is intended primarily for 
Earth orbital trajectories^ it may be used for other types of missions. 
However » a number of options that were included in the original SPACE-A 
trajectory generation program are not presently operational. These op- 
tions were a capability for modeling the perturbations due to thrust, 
a model for lunar libration^ simple models for lunar gravitational ob- 
lateness perturbations, as well as models for the drag of Mars and 
Venus, and a capability for computing on-board observations from a 
vehicle. Although these options are not presently available in the 
modified SPACE-A program, with some effort they could also be included. 



SECTION II 
THEORY 

GENERAL DESCRIPTION OF TRAJECTORY COMPUTATION 

Orbit prediction or trajectory computation is the process of cal- 
culating the position and velocity of a vehicle at any time subsequent 
to some initial time, provided the initial position and velocity of 
the spacecraft are given. 

To accomplish this prediction, one uses the laws of celestial 
mechanics embodied in the differential equations of motion. The 
forcing functions for the equations of motion are obtained from a dy- 
namic model which accounts for the accelerations on the vehicle. A 
reference frame is erected to express the components of the various 
vector quantities, and the equations of motion are numerically inte- 
grated subject to the given initial conditions. Once the vehicle tra- 
jectory is determined, the program can then generate observations. 

The coordinate system used in this program is based upon the 
position of the mean Earth's equator and equinox at 0" 1 January of 
the year subsequent to the initial input time of the program. Coordi- 
nate directions of this frame are inertial with respect to the fixed 
stars; the center of origin of the system, however, may be transferred 
from one central body to another so that the spacecraft motion is 
specified relative to a point mass which itself has a proper motion. 
This reference frame is called the Base Date System or simply the in- 
ertial system. 

Observations made from the Earth are necessarily in a system 
different from the Base Date System. To accomplish transfer between 
various coordinate systems, transformations are provided (see COORDI- 
NATE SYSTEMS AND TRANSFORMATIONS). These transformations include the 
dynamical effects of precession, nutation, and daily rotation of the 
Earth. 



All accelerations acting on the vehicle are specified in the Base 
Date System, The gravitational attractions of bodies in the solar sys- 
tem are functions only of position with respect to the vehicle; conse- 
quently^ the program employs an ephemeris giving planetary coordinates 
relative to the Sun^ and lunar coordinates relative to the Earth, all 
in a Base Date System. A Base Date System is specified for overlapping 
two-year blocks of data, the date corresponding to the middle of the 
two-year file. Specifying an initial time causes the program to choose 
an ephemeris file having as its Base Date the beginning of the year 
following the initial input time of the program. In this way, at least 
one full year of ephemeris information is available before a change of 
reference system is necessary. 

Another acceleration specified in the Base Date System without 
transformation is that arising from solar radiation pressure. Since 
this acceleration is a function of relative position between the Sun 
and the vehicle, its direction is given in the proper frame by using 
information from the ephemeris. 

Other accelerations, such as perturbation* due to Earth oblate- 
ness effects must be transformed through nutation and precession to the 
proper frame. All positions, velocities and accelerations are computed 
in canonical units (i.e.. Earth radii (ER) , Earth radii per hour 
(ER/hr), Earth radii per hour per hour (ER/hr^), respectively) and ap- 
propriate constants are used to transform to other units of measure. 

Vehicle motion is always computed relative to some reference body; 
a planet; the Moon; the Sun. Consequently, the equations of motion 
contain a term which accounts for the acceleration of the reference 
body on the vehicle. The remaining accelerations are usually, but not 
always, much smaller than this primary acceleration and are therefore 
called perturbations. In most cases, they can be regarded as giving 
rise to small disturbances in the orbit determined by the reference 
body and the initial conditions. 



Reference bodies are changed during a trajectory calculation when 
the vehicle leaves the "region of Influence" associated with a partic- 
ular body. Regions of Influence are computed for a body with respect 
to the object of which it is a satellite. Hence, each planet has a 
region of Influence defined relative to the Sun, and the Moon hns a 
similar region defined relative to the Earth, In transferring into or 
out of such a region, velocity as well as position with respect to the 
new reference body must be calculated. 

Since no analytic solution exists for the equations of motion, 
numerical methods are employed to compute the components of position 
and velocity. In the program, a choice may be made between using 
straightforward integration and using Encke's method. The former 
technique, called Cowell's method, is conceptually simple, but suffers 
from precision and machine running time problems. Encke's method, 
although somewhat more complicated, gives dividends in both precision 
and machine efficiency. In this procedure, the Keplerian orbit 
arising from the reference body central force field is taken as a 
nominal trajectory. The perturbation accelerations are integrated 
and the resulting position and velocity increments are added to the 
Keplerian solutions. Naturally, Encke's method is most effective when 
the perturbations are small. 

The present program has the capability to compute a number of 
ground-based observationc • Corrections are provided for the refraction 
of an electromagnetic signal by the troposphere and by the ionosphere. 
Adjustments are computed for errors in elevation angle, range, and 
radial velocity; other angular corrections are calculated from the 
adjustment in elevation angle. 

Description of the basic equations used by the program for the 
dynamic model now follows. It includes derivations of the accelera- 
tions due to the planets. Earth oblateness, drag, and direct solar 
radiation pressure. It also discusses the coordinate systems, numeri- 
cal integration methods, and observation calculations including the 



refraction model used by the program^ 

DYNAMIC MODEL 

The equations of motion for a space vehicle are second-order 
differential equations which describe the accelerations arising from 
the forces acting on the vehicle. These forces can be classified as 
follows: 

(1) Gravitational acceleration of the reference body - primary, 

(2) Gravitational perturbations due to other bodies (e.g,, 
planets) , 

(3) Gravitational perturbations due to reference body 
oblateness, 

(4) Perturbations due to drag, 

(5) Perturbations due to direct solar radiation pressure 
on the space vehicle. 

For orbit determination the reference body primary gravitational 
force almost always predominates the perturbation forces given above. 
The only exception occurs during reentry of a vehicle into the atmos- 
phere, when drag forces may be larger than the primary force of gra- 
vity. Another force which may exceed the primary force of gravity is 
that of vehicle thrust. Since this program is designed primarily for 
orbital computations, this force is not discussed here although it 
can be readily included if desired. 

The simplest gravitational force field is that due to a single 
point mass or equivalently that due to a homogeneous ponderable body 
which is perfectly spherical. In this case the equations of motion 
of a vehicle with respect to the ponderable reference body are: 



R3 



where 

y - G(M-hn) = GM, since m << M 

G is the universal gravitational constant, 

M is the mass of the reference body, 

m is the mass of vehicle, 
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(1) 



R is the vector position of the vehicle with respect to 
the reference body center. 

With Initial conditions Rq and Rq Equation (1) defines a **two- 

body" or Keplerlan orbit, which arises from the primary reference 

body central force field and which may be expressed in closed form in 

terms of its true or eccentric anomaly. 

The above "two-body** dynamic model may be extended to include 

perturbational accelerations acting on the vehicle, in which case 

Equation (1) becomes : 



_ ^+ Pj + P2 + P3 + p^ . (2) 



r3 



where 



?l is the summation of all the accelerations arising from 

planetary, lunar and solar attraction, 
-♦^ 
P2 is an acceleration arising from the oblateness of the 

reference body, 
-♦^ 
P3 is an acceleration arising from drag as derived from 

an appropriate model to be described later, and 

Pi+ is an acceleration due to direct solar radiation upon 

the vehicle neglecting reflected sunlight from the 

reference bodies or planets. 

dt •*■ -^ -► 

The perturbational accelerations Pj, P2, P3, and Pi+ are obtained 

from specific dynamical models whose details are described below. 

PERTURBATIONS DUE TO PLANETARY ATTRACTIONS 

The general expression for the perturbational acceleration, Pj, 
of a space vehicle due to the gravitational influence of the Sun, Moon, 
and other planets (excluding the reference body) is given by: 



J 



(3) 



where 

yj - GMj 

G is the gravitational constant, 
Mj is the mass of the jth body. 



Rvj is the position of the vehicle with respect to the 
Jth body, (Figure 1) 

Rrj is the position of the reference body with respect 
to the jth body I (Figure 1) 

LK is the position of the vehicle with respect to the 
reference body (Figure 1) • 




vehicle 
AR 



jt^ body 



reference body 
Figure 1. Vectors for Planetary Attractions 

If the two terms of the bracketed expression in Equation (3) are 
very nearly equals a loss of accuracy due to round-off errors by 
machine computations will occur. Equation (3) can be written in a 
more convenient form due to Battin ^ -I thereby eliminating this pro- 
blem. Using Equation (3) ^ and Equation (A) below and rearranging terms 
one obtains Equation (5). 

LR - Ryj - Rrj W 

This latter form is used to achieve more computational accuracy. The 
actual programmed equations^ which can readily be shown equivalent by 
substitution are: 

^1 " ^^ f^'^J ^^^"^^ - ^^^ ^^^ 

where 

U[3 + U (3+U)] 



f(U) 



1 ^ (l^U)"^ ^'^ 



m 2 (Rrj + m ' ^ 



W 
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(8) 



The forces Introduced by solar, lunar and planetary attractions 
in orbit prediction about the Earth are usually smaller by a factor 
of 10"^ than that of the primary attraction of the Earth's gravita- 
tional pull. Planetary perturbations (including solar and lunar per- 
turbations) do have a significant effect on a vehicle trajectory over 
long time periods or for deep space probea* 

PERTURBATIONS DUE TO EARTH OBLATENESS 

The fact that the Earth is not a perfect sphere of uniform den- 
sity gives rise to perturbational accelerations due to Earth oblate- 
ness. Formulas are derived below which the program uses to include 
these perturbations in the prediction of near-Earth trajectories. In 
a coordinate system attached rigidly to the Earth, the Earth oblate- 
ness perturbation may be treated as the acceleration of a conservative 
force derivable from a potential. Thus: 

P2 - - VU (9) 

The potential function is given by Equation (12) and the final form of 
the acceleration P2 is given by Equations (26)-(29) and Equation 
(32) below. 

The Geopotential Function 

The geopotential function can be derived from the basic property 
that the potential, U, satisfies Laplace's Equation'" ' ^^' * -I • 

V2u - 0. (10) 

In spherical coordinates. Equation (10) becomes 



V^U 



r2 cos 



3 



- (r2 cos * — J + -^ [cos * -J^J + — 17^^ —J 



(11) 



Although the solution of Equation (11) is not elementary, the equa- 
tion may be solved by applying the method of separation of variables - 
and using Legendre polynomials^ * '^^ -^ . The solution of Equation 
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(11) is given by: 

U - - -^ I (-f)" I {P° (ain ♦) [Cn,m cos tnX + S^^^ sin tnX]} 
n-0 m«o 

where 

y is the Earth's gravitational parameter » 

Re is the Earth's mean equatorial radius » 

r»X»({) are spherical coordinates (see Figure 2)^ 

Pn is the associated Legendre function (spherical harmonics) 
of the first kind of degree n and order m^ and 

Cn^m and 

^n m fi^® numerical coefficients. 

The Legendre polynomials » Py^(x)> are defined by 



(12) 



Pn(x) --T^-^ [(x2 . 1)"]. 



2 n! dx 



.m. 



and the associated Legendre function Pn^^) by 

?^(x) . (1 - x^)^.^Pn(x). 

dx^ 

The spherical coordinates r» A, <{) are defined with respect to 
the fundamental planes determined by the true equator, and the Green- 
wich meridian (see Figure 2). 



Greenwich meridian 



vehicle 




true equator 



Figure 2. Spherical Coordinates 
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The fundamental term in the expression for U is given hy 
m - n - 0. The terms in which m ■ are called zonal harmonics. 
Inspection of Equation (12) shows that these terms vary only with 
latitude, and hence reflect deviations of the Earth's potential from 
a sphere of uniform density that are symmetric around the spin axis 
(e.g., a pear shaped Earth potential can be modeled with zonal harmon- 
ics). Evaluating the 2,0 term for U we have 

" r ^T^^ ^2 (sin ^) C2,o O^) 

and since 

P2(x) - I (3x2 . 1) 

Equation (13) becomes 

" r ^T"^' 2 ^^ ^^"' * - ^^ ^2.0- 
The zeros of Equation (13) are at * 35?25, 




Figure 3. Zonal Harmonic for n = 2, m = 

The sectorial harmonics arise when m - n. The zeros of the 
sectorial terms lie along lines of longitude. The remaining combi- 
nations are tesseral or square in that these terms vanish both along 
a number (m - n) of parallels of latitude and a number (2m) of meri- 
dians of longitude* Figure A provides an illustration of zonal, 
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sectorial, and tesseral harmonic variations for a sample set. 






Zonal 
Harmonics 

n - A 
m - 



Sectorial 
Harmonics 



Tesseral 
Harmonics 



n 

m 



n 
m 



Figure 4. Variations 



In summary, Equation (12) describes the Earth's gravitational 
potential at any point in space; the negative gradient of this poten- 
tial gives the corresponding acceleration. Note should be made, how- 
ever, that the acceleration obtained is in the Greenwich coordinate 
system, and hence must be transformed into the inertial system (tlie 
Base Date System) to be employed in trajectory computation. 

Computation of Acceleration.Duft -ttritBrth Oblateness 

The components of gravitational acceleration are now expressed 
as the sum of terms which have the same general form for any m, n 
combination. First, it is necessary to develop expansions for cos mX 
and sin mX in terms of cos X and sin X, From DeMoivre's theorem, 

cos mX + i sin mX - (cos X + i sin X)™ 

Now, expanding by the binomial theorem 

m _, , 

cos mX + i sin mX - ^ (, )(cos X) (i sin X) , 



k-o 



For m an even number 
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2 
COS mX + i sin mX - I ("j^,) (-l)^(cos X)'""^^(sin X)' 

k-o 



2 ■■• 



k-0 



and for m an odd number 

m-l 
2 

cos mX + 1 sin inX - J] (^j^) ("D^Ccos A)'^"^^(sin A)^'' 

k-o 

k-o 

Equating real and imaginary parts above ^ we have 

M 
cos mX - j; (-1)^ (2'k^ (cos X)"^'^^ (sin X)^^ fiO 

k-o 

sin mX - Y (-1)^(2^^^) (cos X)'"-^^-!^^^^ ;^)2k+] ^^.,, 



k-o 

i_ ^^'^ m! 
where I.J" ™, . ■■ ; ■■ , ' 
^Z-' II (m-Z) ! 



'^^ if m is even 

M [ 

^ if m is oddi 



^^ - 1 if m is even 

and M' - X . 

^ if m is odd. 



From the definition of spherical coordinates: 
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sin 4> 



cos ^ cos X « — 
r 



cos 6 sin X ■ ^ . 



(16) 



where x^ y^ and z are the coordinates of the vehicle in the Green- 
wich coordinate system. 

Using Equation (16) in Equations (14) and (15), we have: 



cos mX ■ 



x 



M 



m 



(r cos (J)) k' 



M' 






'2k^ 



sin mX ■■ 



I (-""(CJw 



2k+l' 



(y) 



(r cos (J)) k-o 
Substitution of Equations (17) and (18) into Equation (12) yields: 

^ ■ - r I Lf J ^ t'^ — ::i ^(''» y^ 



n"0 



m-o (^ cos ({)) 



where 



M 



G(x, y) -Cn.ni I ("D^l^l^ (x)°^-^^y) 
k«0 

. c V / T\^r ro ^/ sni-2k-l. v2k+l 



k«0 



From the definition of the associated Legendre function, we have 

m 



rjm 



(T) - ^^ - ^ ^ ^-=-[(12 - 1)"]. 



2" n! 



dT 



m+n 



Setting T « sin <fr and dividing Equation (19) by cos 0; 
^n («^" ») . _JL_ J^ [(t2 - 1)"] 



cos™ '^ 
n 



on , J in+n 
2 nl dx 



Setting U - J] I Un^nt ^^'^ using Equation (20) we have 
n»0 n"0 



(17. 



(18) 
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(]9) 



(20) 
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dT 
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- U R. 



m+n+l ^n , 
r 2 n! 



•^ (t2 - 1)" C(x, y) 



dT 



tn+Ti 



The gravitational acceleration is computed by taking the negative 
gradient of the potential function. The general derivative in the 
gradient will be taken with respect to C, where ? takes on the 
values of x , y ^ and z • 

Defining 



Am,n 



3U 



n>in 



n 



and carrying out the indicated operations using 



H r 



and 



il _ J- fil i 1 riS. _ 2l 



yields: 
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dT 
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r L3C r^J , m+n+l 



dT 



d"^:^ (,2 . ij 



dT 



m+n 



n 3G(x, y)" ! 
3C J 



where 
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3C 



k-o 
in-2k 2k- 1 ayi 



3C 



+ 2k X- - y~" -||}] 



3C 



M' 

r r , ,vk /■ m If , ~, ,N ra-2k-2 2k+l 9x 
+ Sn.in I [(-1) (2k+iJi(ni-2k-l) X y — 

k-0 
+ (2k+l) x' 



in-2k-l 2k 9yi, 



(21) 



(22) 
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It can be ihown by applying the binomial theoran and induction 
that ^ [n - 4] 

,,i ^ ^ / ^ ^^ ^k^ (2n-2k-i)! ^^^ 

where [n - *] indicatea the integral part of n - T« Hence, 

[221] 

^ m+n ^ ^ i ^ ^^ IjcJ (n-2k-m)! ^^^ ^"^ 

dT k-O 

and , 

[2=£=i] 

jPi-t-n+l n V / 1^k fni (2n-2k'> ! . .n-2k-m-l ,,,, 
dT k«o 

where t ■ sin <> ■ * » Therefore, Equation (21) may be rewritten as 



y R 



rEIILi 

n L— JT-J 



a5 *" ^ ^% fj f / l^^ (^\ (2n-2k) ! |^z^n-2k-ml 

2 n! r k«0 

(iasiii 5 OCX. y) - i£i|iii) - 

k-0 

X {2J^ [|| - i^] }] (25) 

y 

By separately considering the case for n-m odd and for n-m even, 

it may be shown that Equation (25) reduces to the following expres- 

sion: ^ ^ l^] 

.5 _ » ^e" r ,. ..k (2n-2k)! .z^n-2k-m-l1 .. .... 

^"•"^ ■ - ,n nH^+1 ^ ^^-^^ k! (n-k) ! (n-2k-m) ! U^ ^' ^^*^^ 

{((2n-2k+l) f .lazuli) c(x. y) " f -i^^l 
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where 



2x _ rl if C - X ^ rl if C - y ^z rl if C - 
3C " ^0 if C »* X* 3C " ^0 if C ?< y* 3C " ^0 if C ?< 



«"d M 



,k w! , .ra-2k, .2k 



-+ 



k ml 



Sn,m I [(-1) 



(2k+l)! (m-2k-l)! 

k-0 



(27) 



G(x. y) - Cn^„ I (-1) 2kl"n.-2k)! (^)""' (y)' 

k-0 

M' 

. o V / i\^ SJ / xni-2k-l .2k+l -^-. 

+ Sn.m 1 (-1) (2k^l)! (Tn-2k-l) ! ^^^ ^^^ ^^^^ 

k-0 

r /o 4r 4 (^ - 1 if m is even 

_ I m/2 if m is even m' ■ / 

l.-^- if in is odd l^T" ^^ ^ ^® ®^^ 

^ — ^n,in L K i; 2kl (m-2k)l ^ ^"^ ^^mx; kv) 3^ 

k-0 

+ 2k(x)'"-2Ny)2'^-l|Z}] 



// ^1 -.N/ Nm-2k-2. v2k+l 3x ^ /oi^i\/ xm-2k-l, .2k 3y 1 , .^-. 

X { (m-2k-l) (x) (y) rr + (2k+l) (x) (y) -i}] (29) 

d^ d^ 

It is important to note that when 2 - 0, a special procedure 
must be used since the term (— )"^ arises in Equation (26) for n-m 
even; and (— j^* when n-m is odd. 

When n-m is even. Equation (26) reduces to one term, when z = 0. 

n-m 

(30) 
2 -/* V" 2 



^ n— m ^ ^ 

C - V ^e If ,. 2 (n-Hn) 1 ^ Kn-hiH-PS r. ^l 

^n.m - ^, ^,^^1 ^(-1) ^^_^^ J 1^ ,2 - acj • 



Similarly, when n-m is odd. Equation (26) may be reduced 
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An,ta - for 5 - x» y 

A^,, - ^^^, ((-1) ' 2 ^ teiill—^ £ (31) 

In summary, the program computes the acceleration of a vehicle 
due to the Earth in the Greenwich coordinate system by the formula: 

® n ^ 

P2 - I I (Ajl^ni ^ + Aj;^ni J + A^^m k) (22) 

n-0 m-0 

"t "^ '^ 
where i» j» k are unit vectors in the Greenwich coordinate syHtem. 

In Equation (32) the range of indices (n,ra) are riestricted by 

SPACE-A to the following limits: 

2 0, 1, 2 

3 0, 1, 2, 3 

4 0, 1, 1, 3, 4 

5 0, 1, 2. 3 

6 0, 1, 2 

7 0, 1 

8 

9 
10 

The fundamental term is given by the (0, 0) combination and 
would be the only terra present if the Earth were of uniform density 
and a perfect sphere. Since the center of the coordinate system is 
taken to be at the center of mass, it can be shown that the (1, 0) 
and (1, 1) combinations vanish*- '^' -^ . 

In Equation (32), A is given by (26) when z 7^ 0, and when 
z = 0, by Equation (30) or Equation (31) depending on whether n-m 
is even or odd. 
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PERTURBATION DUE TO DRAG 

Accurate simulation of an artificial satellite or other space 
vehicle trajectories requires consideration of vehicle deceleration 
resulting from atmospheric drag. A number of planets (e»g.p Earth, 
Venus, Mars and Jupiter) have sufficiently dense atmospheres to re- 
tard the motion of a vehicle within their atmospheres; however, this 
discussion confines itself to the atmospheric model of the Earth. 

An analysis of drag must take into account the particular mis- 
sion of the vehicle, e.g., low eccentricity, orbit, reentry, or fly- 
by, since vehicle mission determines what portion of the atmosphere 
it is necessary to include in the model. 

The following discussion describes the general equations used 
for drag computation, some of the problems involved in simulating the 
Earth's atmosphere, and the effects of certain simplifying assumptions. 

Drag Equations 

The program uses two different formulas for computing the vehi- 
cle deceleration, P3, due to drag. For a relatively dense atmos- 
phere where the assumption of continuum flow is valid the following 

well known equation is used: 

-V 1 CnS H. 

P3 - - Ci [^ P -jjp] Va Va (33) 

where 

p is the atmospheric density at the vehicle, 

CD is the drag coefficient of the vehicle (dimensionless) , 

S is the effective surface area of the vehicle, 

m is the mass of the vehicle, 

Va is the vector of the velocity of the vehicle with 
respect to the local atmosphere, 

-► 

Va is the magnitude of Va» and 

Ci is a constant used to convert the above expression 
into the units used by the program. 

Equation (33) is used to compute the acceleration of drag in the 
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lower atmosphere from kilometers up to 120 kilometers, although the 
model for the lower atmosphere from which the values of p and C^ 
are obtained may be extended up to 210 kilometers with some loss of 
accuracy. The basic model used in the lower atmosphere is the U, S. 
Standard Atmosphere 1962^ -' ' and its details are discussed belov;. 

As the atmosphere becomes more diffuse, the mean free path length 
(average distance between impacts of air molecules) increases. At 
110 kilometers, mean free path length is roughly one meter and at 130 
kilometers may become as large as ten meters -* - When the maan free 
path length exceeds the dimensions of the vehicle, the assumption of 
continuum air flow is no longer applicable. In such a diffuse atmos- 
phere where all collisions are essentially two-body collisions, the 
air flow is referred to as free molecular flow. 

Ketchum -' has derived the following formula for the magnitude 
of drag deceleration in free molecular flow: 

2R>ir ^ S 



IP3I -t(1+^][p Cavr] Va (34) 



8 

where 

R the radius of the vehicle, 

X the mean free path, 
Cav the average velocity of particles in the medium. 

Ketchum is uncertain as to the validity of the (1 + 2R/X) term 
in Equation (34). In the high atmospheric region the assumption that 
X >> R is usually justified, except perhaps below 140 kilometers. 
Therefore, the program actually uses Equation (35) in computation of 
the drag in the high atmosphere. 

o m 
where, 

S is the effective surface area of the vehicle, 

m is the mass of the vehicle, 

p is the atmospheric density at the vehicle, 

Cav is the average velocity of particles in the medium, 
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- C2 [t P Cav r] Va (35) 



Vg is the vector of the velocity of the vchic]c with 
respect to the surrounding atniosphere, and 

Cp is a constant used to convert the expression into 
the units used by the program. 

Equation (35) is used to compute the deceleration of drag in the 
upper atmosphere from 100 kilometers up to 1^050 kilometers, although 
the model for the upper atmosphere may be extended up to 2>000 kilo- 
meters* The basic model used for the upper atmosphere is that due to 

r9i 

Harris & Priester'" -^ . The values of p and Cav used in Equation 
(35) arc derived from this model, the details of which are also rlis- 
cussed below. 

Notice in both Equations (33) and (35) that the direction of the 
drag force is in a direction opposite to the velocity with respect to 
air. In addition, both formulations assume zero lift and assume that 
the angle of attack of the vehicle is zero, i.e., that the vehicle 
velocity relative to the air mass is in line with the vehicle lonei- 
tudinal axis. 

It is readily seen that the formula for drag in the upper atmos- 
phere. Equation (35), differs from the equation of drag in the lower 
atmosphere, Equation (33). Furthermore » in the region between 100 
kilometers and 120 kilometers there are disagreements between the two 
models. For example, the value of density predicted by the low atmos- 
pheric model (U. S. Standard Atmosphere 1962 ) nt ]2n kllorcters is 
34% higher than that predicted by the high atnospheric rrdcl (Varrlr.- 
Pricstcr) at the same height. 

The present progrnn achieves a compromise between tncpe tvo 
models by treating the region from 100 to 120 kiloretcrs ns a trr*^?-*- 
tion region in which a weighted average is taken between the drar 
values computed by the two r.cthods and gradual]y slides the weight 
fron unity for free molecular flow and zero for continuum flow at 120 
kilometers to unity for continuum flow and zero for free rolecu]ar 
flow at 100 !;ilometers. 
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Because of the uncertaintiei in the atmospheric models and be- 
cause of the approximation! made in the analysis the computation of 
drag deceleration is probably accurate to + 5% in the lower atmospi.ert 
and is less accurate in the upper atmoiphere. 

Variables Used in Drag Computation 

Vehicle Mass 

In the most general case^ the vehicle mass in the drag equations 
should be considered as variable with time. In the orbiting case or 
the fly-by case, a step change in mass representing the separation of 
a landing craft is conceivable, A long-term steady-state mass flow 
rate, however, would probably be small. 

For the reentry case, if the reentry vehicle is of the heat-sinK 
type, the mass would be constant. For an ablative nose cone (i.e., 
one which loses mass when moving at high speeds due to friction) the 
mass flow rate is a function of the drag. For ballistic missil*^ ap- 
plications, this mass change is usually ignored. In any event, such 
changes in mass represent a small error in the location of the impact 
point. Therefore, the program treats the mass as an input constant. 

Surface Area 

The effective surface area term S in the drag equation is not 
simply the cross-sectional area of the vehicle. The vehicle, in 
passing through the air, produces a shock wave which skirts the mis- 
sile thereby placing the effective cross-sectional area at a point 
somewhat close to the nose. Since the shock wave changes with air 
speed, so does the effective cross-sectional area. In practice, S 
is made constant and any variation with speed is included in the co- 
efficient of drag. 

Velocity with Respect to Air 

The velocity of the vehicle is available in an inertial coordi- 
nate system (Base Date System), 

The vehicle velocity with reapect to the moving air mass, Va, 
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in the same coordinate systemi is obtained by subtracting the velocity 
of the air mass from the vehicle velocity, A good first approximation 
to the velocity of the air mass is obtained by assuming the air mass 
to be rigidly attached to the rotating planet. 
From these considerations we have: 

Va - R - «' X R (36) 

where 

Vq is the velocity with respect to surrounding atmosphere ^ 

R is the position vector of the vehicle in the inertial 

frame t 
-► 
R is the velocity vector of the vehicle in the inertial 

frame 9 and 

J] is the vector of angular rotation of Earth expressed 
in the inertial frame, 

A better approximation could be obtained by including the effects 
of wind velocity. The purely local effects have to be neglected , but 
the long-term horizontal effects are known as a function both of posi- 
tion on the Earth's surface and of altitude. The effects of the wind 
velocity's direction (independent of altitude but dependent on lati- 
tude and longitude) and magnitude (strongly dependent on altitude , 
less strongly on latitude ^ and least on longitude) would have to be 
included. The error made by neglecting Earth winds is about 1,500 
feet at impact for a typical ICBM mission. It should be noted that 
winds are of importance only in the Earth's lower atmosphere, mainly 
for the reentry case. In the present program the effect of winds is 
neglected. 

Drag Coefficient 

The drag coefficient, Cd, is sometimes considered to be a con- 

CdS 
stant. Very often the ballistic coefficient, B - — - , is used in 

" m 

analysis in place of the drag coefficient, mass, and effective sur- 
face area, A much more accurate representation for Cd is obtained 
by considering it to be a function of Mach number, Mach number is 
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defined to be: ^ 

V. 

(37) 







M - 


c 


e, 








M 


Is the Mach number. 







\Vg^\ is the speed with respect to the surrounding alr^ and 

c is the speed of sound in the surrounding air* 

The speed of sound is a function of altitude obtained from the 
low atmospheric model (see below). It should bs noted that as altitude 
increases, the atmosphere becomes rarified to the point that the speed 
of sound loses its physical significance. 

In the program Cd la tabulated for about AO different Mach 
numbers. These numbers are denser for speeds below Mach 2 than those 
above and are very dense in the region around Mach 1, For intermedi- 
ate values of Mach number, linear interpolation is used. 

Inadequate knowledge of the drag coefficient is one of the major 
sources of inaccuracy in the simulation of drag. Since drag coeffi- 
cient is a function of Mach number, drag coefficient data has been 
obtained by wind tunnel measurements made at a range of Mach numbers. 

Air Density 



In the region below 120 kilometers the air density, p» is a 
function of altitude obtained from the low atmospheric model and is 
computed from a stored table (see below). 

In the high atmosphere p is obtained from the upper atmospheric 
model of Harris-Priester and is considered to be a function of alti- 
tude, local solar time, solar flux, and latitude. Tables are provided 
in the program for its computation (see below). 

Mean Particle Velocity 

The mean particle velocity, Cav* i* of importance only when 
the vehicle is in the upper atmosphere where the assumption of free 
molecular flow is valid* From Equations I«3«4-l and 1,2,6-1 of 
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Reference L7], C is given by 

where 

T is the absolute temperature of particles (*K), 

m is the mean molecular weight of medium (gm) , 

k is a constant of proportionality (.IAS). 
The values of T and m are given in the Harris-Priester model of 
the upper atmosphere* Cav can, therefore, be obtained directly from 
stored tables. 

Lower Atmospheric Model 

Drag in the lower atmosphere (below 120 kilometers) can be large 
and a vehicle entering this region will usually be slowed down suffi- 
ciently to be quickly captured by the Earth. Thus, the lower atmos- 
phere is primarily of concern in the reentry case. 

Data for an average model have been well established for the 
lower atmosphere. There are five sources for these data: U. S< 
Standard Atmosphere 1962 ; COSPAR International Reference Atmosphere 
fCIRA), 1961; COESA Table for Tropical Latitudes, 1962; ARDC Model 
Atmosphere 1956, 1959. Table I shows the density deviation (in 
percent), as a function of altitude, of each of the others from the 
U, S. Standard Atmosphere values. From the table, it is evident that, 
except for the COESA tables, there is good agreement between the vari- 
ous tables at low altitudes. Note that the U. S. Standard Atmosphere 
and CIRA tables are in excellent agreement all the way to 120 km. 
(400,000 feet). 

The lower atmosphere is characterized by seasonal, diurnal, and 
latitude variations; however, none of these is sufficiently well 
documented. The only effect of omitting them is that the impact point 
of a re-entering body would be slightly different. It was estimated 
in 1958 that the standard deviation for a heat-sink type nose cone 
used in the ICBM application is only about 0«5 nm, 
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Table I 



Comparison of Sources of Density Data 
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Three tables with 32 values each are stored in the program for 
the lower atmospheric model. The first contains 32 values of alti- 
tudes in kilometers, the second contains 32 values of the logarithm 
(base 10) of the density, p, in gm/km^, and the third contains 32 
values of the speed of sound, c, in km/sec. The n^h entry in the 
altitude table corresponds to the n^h entry in the table for logjoP 
and c. Intermediate values are obtained by linear interpolation, 

Tliese three tables are obtained from the U, S, Standard Atmos- 
phere 1962 > in which densities and speed of sound at altitudes below 
90 km are listed. Above 90 km the speed of sound was calculated from 
temperature and mean molecular weight data which are directly avail- 
able. 
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The value of altitude (in Earth radii) above an ellipsoidal 
Earth is obtained using a formula found in Baker^ -^ . 

h - r - 1 + f sin^ 4, + -L- (i « i) gin^ 2(J) (39) 

where 

h is the height of the vehicle (ER) , 

r is the distance of the vehicle from the geocenter (ER) , 

(^ is the geocentric latitude of the vehicle, 

f is the flattening constant of the Earth ■ -,p . . 

For programming purposes Equation (39) is put into the more convenient 
^o^^* o 1 2 -i. v2 

, , ,, , liiiiiiKMoili 

^ •• ^ - ^ + ^ 298.3 R2 (40) 

where , 

h is the altitude of the vehicle (ER) , 

X|y»z are the position coordinates of the vehicle in the 
Greenwich coordinate system (ER), 



/ 



9 . o . o is the distance of the vehicle from 

X 1 V T z 

^ the geocenter (ER) • 



The position coordinates of the vehicle are available from the pro- 
gram and h is multiplied by 6378.165 to convert its units into 
kilometers. 

After h is obtained, p (p - 10 ^^10^) and c are obtained 
from the low atmosphere tables. Then the air velocity, Va and the 
drag coefficient, Cp, are computed as described above. Finally, 
P3 is computed according to Equation (33), 

Upper Atmospheric Model 

Models of the Earth's upper atmosphere (above 100 kilometers) 
must take into account solar activity. There is evidence that solar 
activity occurs cyclically at periods of 27 days, 6 months, 1 year 
and 11 years. 
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Theoretical modals do not exist for the 27--dayi 6-month, and 1- 
year cycles. Diurnal variations, If any, of the models for these 
cycles are not known* Investigation of the 11-year cycle (corres- 
ponding to the sunspot period) in solar flux led to the Harris- 

I 9 111 
Priester model of the upper atmosphere. This model*" * -^ has 

diurnal and solar flux variations. 

r9i 

The Harrls-Prlester model"" -■ lists the density, absolute tem- 
perature, and mean molecular weight of the atmosphere as a function 
of altitude, solar flux and local solar time. The mean particle 
velocity can be found by use of Equation (38). Note that at the 
North and South Poles local solar time Is undefined. 

The upper atmosphere has a delaying effect on solar radiation. 
It takes several hours for the Sun's heat to pass through the atmos- 
phere and reach the Earth's surface. The Harris-Priester model is 
based on densities computed at the Earth's equator. Intuitively, it 
is expected that it will take longer for the solar flux to reach the 
poles as opposed to the equator. Therefore, it is considered that 
there is an effective variation of solar flux with latitude. This 
variation is implemented In the program by applying the Harris-Priester 
model at the equator and a stored table of "twilight" densities at 
the poles . The cosine of the latitude of the vehicle is used as a 
weighting factor to Interpolate between the two sets of data. 

The Harris-Priester upper atmospheric model has been incorporated 
in the program by means of a table look-up procedure.. Two tables are 
used in the program. 

The first table lists the logarithm to the base 10 of the density 
in gm/km^ times the mean particle velocity in km/sec (logjoP ^av) at 
the equator as functions of altitude (km), solar flux (watts/m^ - Hz 
at a wavelength of 10.7 cm) and local solar time (hrs). These values 
have been tabulated for 16 values of altitude, 4 values of solar flux 
and 13 values of local solar tine. Hence, this table has 832 entries. 
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The second table lists logjo pCav in the same units at "twi- 
light*' for the modeling of the polar region as a function of altitude 
and solar flux. This table has 16 entries for altitude^ A entries 
for solar flux, or 64 entries. 

For intermediate values of the input variables (altitude, solar 
flux, local solar time) a linear interpolation is used to obtain the 
output, logio pCav* This method gives fairly accurate results since 
oCav is nearly exponential. The value of solar flux is determined by 
input data. The value of altitude is computed according to Equation 
(40). The value of local solar time is computed from the x and y 
coordinates of the vehicle and the Sun in the inertial (Base Date) 
coordinate system (see Figure 3). 
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Figure 5. View of x-y Plane from Above 
Thus we have: 

local solar time-{(e2 - Oi + tt) -122-+ 720'*}(mod 360**) 



where 



ys 



6i ■ tan"^ f"^l where ir < 0i < - tt 



Xs 



Go ■ tan"^ ("•] where tt < e^ < - tt 
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Xy, yv are the vehicle coordinates in the inertial system, 
xs» ys are the Sun's coordinates in the inertial system. 

Local solar time given above is then divided by 15 to convert Its units 

to hours. 

After a value of .logjo pCav baa been obtained from both the 

equatorial table and the "twilight** polar table an interpolation is 

done to obtain the final value of pCav using the cosine of latitude 

as a weighting function* Hence: 

L " Lp + cos ^ (Le - Lp) (^"i 



/y + v^ 

"/x2 + y2 + 



cos ^ - /.J ^ y2^ , 9 (^3) 

pCav - 10^ (4/.) 

where 

L is the final value of logjo PCav» 

Lp is the value of log^o P^av obtained 
from the twilight table, 

Le is the value of logjg PCav obtained 
from the equatorial table, 

x,y,z are the coordinates of the vehicle in the inertial 
(Base Date) system and ^ is the latitude of the 
vehicle. 

Once pCav has been found, the vehicle area and mass (S and m, 

standard program inputs) as well as air velocity V^ (see Equation 

(36)), are obtained and the value of drag deceleration is finally 

computed according to Equation (35). 

PERTURBATION DUE TO DIRECT SOLAR RADIATION 

Solar radiation exerts a pressure on the intercepting surface of 
a vehicle. Orbiting planetary vehicles, having a large area to mass 
ratio are subject to noticeable perturbations due to solar radiation 
pressure. In fact, for orbits above 500 miles the solar radiation 
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[121 
perturbation is usually more significant than that due to drag^ -'. 

The vehicle acceleration due to solar radiation pressure depends on 
the area to mass ratio of the vehicle as well as the intensity of the 
Sun's incident power at the vehicle and the fraction of solar illu- 
mination on the vehicle* The illumination factor is considered in 
three distinct regions in the following analysis: full sunlight^ 
penumbral illumination^ and no illumination (i.e.^ umbral region). 

In computing the solar radiation perturbation, Pi^, this analy- 
sis neglects the dispursive effects of planetary atmospheres which 
complicate the geometry of the umbra and penumbra. The analysis also 
neglects the effects of reflected sunlight from the reference body 
or any other planet. 

Acceleration Due to Solar Radiatj^ y^ yyggpni^e 

An expression for the acceleration due to solar radiation pres- 

[121 
sure given in Wolverton ■' is: 



rAK-Sl-l Iss^ 



where J 

p is the illumination factor, 

q is the reflectivity coefficient, 

A is the area of vehicle pertaining to solar 
radiation pressure, 

m is the mass of vehicle, 

Lo*3.86 X 10^^ watts, the total power output 
of the Sun (+ 3%) , 

c is the speed of light, 
Rsv is the vector from the Sun to the vehicle. 

The SPACE program assumes a reflectivity coefficient equal to 
unity. This assumption may be altered by changing the area. A, by 
an appropriate factor. Simplifying Equation (45) and putting the 
above variables into units used by the program, one obtains: 

^ - p Ir) cp fe («> 



^m' * *^sv 
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where > 

p Is the illumination factor, < p ^ 1, 
A is the area pertaining to radiation pressure (ft*^), 
n is the mass pertaining to radiation pressure (lb-masses), 
Rsv Is the vector from the Sun to the vehicle (ER) , 
Cp "1.04819 X 10-^ (liRi-fti^ ^^^ ^^ ^ constant used to convert 

the expression into units used by the program and to ac- 

LO 
count for the 7— factor of Equation (45). 
^ffc 

Equation (46) Is the basic equation used by the program for com- 
puting the solar radiation perturbation acceleration. 

The Illumination factor, p, is obtained from the relative geo- 
metry of the vehicle with respect to the Sun, the Hocn and the planet-;: 
p - 1 in full sunlight, 
p - in the umbral region, 

< p < 1 In the penumbral region. 

It Is possible that a vehicle may lie within the penumbral region 
of two bodies at the same time, e.g., the Earth and the MooPv In 
this case, only the penumbral Illumination factor due to the reference 
body is computed* Errors introduced by this assumption are extremely 
small: first, because for most vehicles of interest the solar radia- 
tion perturbation is Itself small (usually less than 10 ^ times the 
acceleration of the reference body except for low density balloons) ; 
second, because only a short time is spent in the penumbral region of 
the reference body by an orbiting vehicle; and third, because the 
incidence of simultaneous penumbral obscuration Is rare. 

Since the geometry used In calculating the illumination factor 
Is the same as that used for eclipse information, the portion of the 
program which computes the solar radiation pressure perturbation also 
computes the times at which the vehicle enters or leaves the umbra or 
penumbra of the celestial bodies. The geometry used In calculating 
the vehicle illumination factor, p, is described below. 
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Vehicle Illumination Factor 

Figure 6 illustrates the geometry of vehicle illumination, neg- 
lecting the effects of atmospheric refraction. By similar triangles 
it is seen that the height of the umbral cone (hu) and penumbral 
cone (hp) are respectively given by: 

hu - ^ - ^^ (A7) 

Rp 

hp . ~^^ (A8) 

^ + 1 
Rp 
where, 

Rsp is the distance from the center of the Sun to the 
center of the reference body, 

Rs is the radius of the Sun, and 

Rp is the radius of the reference body« 

Next, criteria are developed to see in which region the vehicle 

lies, i»e», full sunlight, umbra, or penumbra. First, consider 

Figures 7 and 8 and the definitions and relations below* 

Rsp is the vector from the center of the Sun 
to the center of the reference body, 

-► 

R is the vector from the center of the 

reference body to the vehicle, 

? - hu |2£ (49) 

Ksp 

-+■ 
•^ , ^SP 

Q " - "^P R^; (50) 

From Figure 7 we get the following: 



cos 



(52) 



\t -?\ Rsp 
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POSITION OF VEHICLE 



SUN- 
LIGHT 




APEX OF UMBRA 



Figure 7. Umbral Region Geometry 



SUNLIGHT 



APEX OF 
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POSITION OF 
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Figure 8. Penumbrol Region Geometry 
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Froji Flgurta 8 we get; 



cos B - /l - [^f (53; 



^hp 



cos B^SLzJLLhz i5L) 

1^ - Q! Rap 

Nov^ if the scalar product R * P is positive , then the vehicle 
lies on the side of the reference body away from the Sun, In this 
ciise, if cos g > and if ex ^ A, the vehicle or satellite lies 
in the unbra and p is set equal to zero. Another test is: 

if I cos a| >, I cos AJ (^"'^ 

ihen p ■ 0| 

i.e., the vehicle lies in the umbra. 

If the vehicle does not lie in the umbral region one can test to 
see if it lies in the penumbra. Thus, If cos 6 > and 6 ^ B the 
vehicle lies in the penumbra, or equivalentlyi 

if I cos 6 1 >^ I cos b| (5^) 

then < p < 1 

j^i.e.^ the vehicle lies in the penumbra). 

If the vehicle does not lie in either the umbral or penumbral 
region, then the reference body does not obscure the Sun's rays, A 
check can then be made to see whether any other celestial body blocks 
or partially blocks the solar radiation; and if not, the illumination 
factor is set equal to one, p ■■ 1, 

Penumbral Illumination Factor 

If the tests of Equations (55) and (56) above indicate that the 
vehicle lies in the penumbral region, then the illumination factor, 
p, is computed according to the formulat 



P 



aex 



e 



s* 



(57; 
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where A£x Is the angular area subtended by the exposed portion of 
the solar cap at the vehicle position^ and 

6s is the total angular area of the solar cap at the 
vehicle position. 

The general approach used by the SPACE program for computing the 
penumbral Illumination factor Is somewhat more detailed than chat 
given In most references* It consists In projecting the solar disc 
(or cap) and the cap of the obscuring planet (or moon) on to a great 
imaginary sphere whose center Is at the vehicle position. The rela- 
tive angular areas of the caps are computed and the Illumination fac- 
tor Is given according to Equation (57), While the theoretical de- 
velopment is somewhat Involved^ It leads to simple closed form alge- 
braic expressions convenient for use on a computer. 

Firsts consider Figure 9 which shows a sphere of radius R^ re- 
presenting the Sun or a planet » and a point P at a distance Jl + ^ 
from the center of the sphere. The apparent angular area of the 
sphere as seen from point P is equivalent to the angular area of 
the sphere's cap projected onto a projection sphere of radius a. 
Therefore the following Is truei 

angular area 6 - 2Tr [1 - cos if] 
- 2Tr 



1 Aa + 2R) 
L ^ + R - 



0p » 2iT 



/h(h + 2Rp) 
1 - 



h + Rp J a 



C^S) 



- 27r ^ . (59) 

a 

where a Is the radius of the projection sphere ^ and 
H is the height of the projected cap. 
The total angular area of a planet (6p) or of the Sun (65) at 
the position of a vehicle can be obtained by substituting the appro- 
priate values of I and R Into Equations (58) and (59) and simpli- 
fying. Thus 9 we have: 



27T - (60) 
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POSITION OF VEHICLE AT 
WHICH THE ANGULAR AREA 
OF THE SOLAR OR PLANE- 
TARY CAP IS MEASURED 



DISC AND CAP OF THE 
SUN OR PLANET 




SUN OR PLANET 



SPHERE OF PROJECTION 



PROJECTED CAP OF SPHERE 

Figure 9. Geometry Used to Measure the Angular Area 
of the Solar or Planetary Cap 
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2n 



- A^ 



2n f (61) 



where 

6p Is the angular area of the planet (or moon) at the vehicle 
position^ 

h Is the height of the vehicle above the planet's surface^ 

Rp Is the radius of the planet ^ 

6s Is the angular area of the Sun at the vehicle position^ 
Rsv Is the distance from the center of the Sun to the vehicle^ 

Rs Is the radius of the Sun^ 
H Is the height of the solar cap» 

H* Is the height of the planetary cap^ 
a Is the radius of the Imaginary sphere of projection. 
Equations (60) and (61) establish the relative size of the solar and 
planetary caps as seen from the vehicle. 

Next^ It Is desired to determine the angular area of the exposed 
solar cap A£x> '^o do thls^ consider Figure 10. In this figure the 
vehicle Is positioned at the center of a great Imaginary projection 
sphere of arbitrary radius ^ a. The relative positions of the Sun 
and the obscuring body or planet are shown as well as the projections 
of the solar and planetary caps onto the great sphere. The equator 
of the great sphere Is constructed to be coplanar with the vehicle^ 
the center of the Sun, and the center of the planet. A great circle 
Is also constructed perpendicular to the equator and passing through 
the Intersection of the caps; this circle will be used to define one 
of the limits of Integration In the computation of the angular area 
of the two lunes which are formed by the great circle. 

For the calculations to follow^ three coordinate systems are de- 
fined and Illustrated In Figures 10 and 11. The coordinate systems 
employed are: 
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-PLANET 



SUN 



SOLAR CAP OF 
TOTAL AREA 65 




PLANETARY 

CAP OF 
TOTAL AREA 

ep 



GREAT CIRCLE 



Figure 10. Geometry of the Intersection of the Solar 
and Planetary Caps 




Figure 11 . Definition of Coordinate System 
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(1) A systttm (x^y^z) for th« solar cap» 

(2) A eyetem (x\y\2') for the planetary cap, 

(3) A ayatetn (x'\y'\z") for the great circle. 

All three coordinate systems ar« orthogonal and have^ In common^ the 
same y axisi The relations bstween the coordinates are giv«n in 
Equations (62) and (63) t 



f *^ 




x" 




z" 


- 


y" 





cos Qq 


sin 


6C 





sin Qq 


cos 


ec 








( 


D 


1 


■in Qq 


cos 


ec 





- cos Qq 


sin 


ec 








( 


3 


1 



(62) 



(63) 



where 



6c is the angle from the center of the solar cap (z-axis) to 
the center of the planetary cap (z*-axis), 

6g is the angle from the center of the solar cap (z-axis) to 
the great circle (x"-axis). 

By inspection of Figure 10 one finds that: 



R . R 



cos @c 



MS. 



^sv 



(6A) 



where 

< ec i •J 

Next find the points of intersection of the planetary and solar 
caps X - Xp, y - yp, and z - zp. Referring to Figures 9 and 10 
one sees that the equation of the small circle of the solar caps 
(having a cap height of H) ist 
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x^ + y2 + z^ - a^ 



z - a - H 

•*. x^ + y2 - H (2a - H) 

Similarly, che equation of the small circle of the planetary cap 
having height H' is: 

(x')2 + (y')2 + (zt)2 . ^2 

z' - a - H' 
x'2 + y»2 . ^f (2a - H') 

Using Equation (62) we can write Equations (67) and (68) in terms of 
the (x,y,z) coordinates. 

z' » X cos 6^ + z sin 6c ■ a - H' 

x'2 + y»2 . (^ j,Qg ec - z sin ec)^ + y - H'(2a - H') 



(65) 
(66) 



(67) 
(o( 



(69) 
(70) 



From Equations (65) and (69) we find: 



(a • 


- H') - (a - H) cos 9c 


X - Xp - 


sin e^. 



Using this result and Equation (65), Equation (70) gives: 



y - yp 



2 H'(2a - H') sin^e^ - [(a - H') cos Bq ' (b - H) ] 2 

sin^e^ 



And, of course, Equation (65) can be rewritten: 



z * Zp * a - H 



Equation (72) requires some interpretation. Figure 10 illustrates a 
case where the solar and planetary caps intersect. It is possible 
that the caps are tangent or that the planetary cap lies within the 



(71) 



(72) 



(73) 
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solar cap. To determine whether an intersection exists, the value 



of yp2 
(1) 
(2) 
(3) 



of Equation (72) may be used as a discriminant. 

If yp^ > 0, two intersections exist. 

If yp^ - 0, the caps are tangent. 

If yp^ < 0, the caps do not intersect and the 

exposed solar angular area is given by 

Aex - 9s - 9p« 

For the great circle passing through the intersection of the 
caps z" - 0, z - zp, X ■ xp. Hence, Equation (63) gives: 

z" - - X cos 0Q + z sin Oq " ^ 

x ^P 

Qi- ^ m — . tan Gp 

Z Zp ^ 

thus from Equations (71) and (73) 



tan 0Q - 



(a - H') - (a - H) cos 6^ 



(a - H) sin 9^ 



In a similar manner we may also find the equation for 9q', 
the angle between the center of the planetary cap and the great 
circle^ 

tan eg' - ^ 

everywhere on the great circle. 

Using Equation (62) evaluated at x ■ xp and z ■ zp we 
finally obtain: 



(a - H') cos e^ -> (a - H) 
"" ^g' - (a - H») sin Be 



In computing Aex» there are five cases of interest and these 
cases are illustrated in Figure 12. If yp^ < one obtains case 
V where Aex " % - 6p« If yp^ * 0, however, the caps intersect; 
and before Aex can be determined, one must calculate A| and 



(74) 
(75) 



(76) 



(77) 



(78) 
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Case I. A£^^ ■ Ag - Ap 



Case II, A£x • ^s * ^s * ^P 
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Case III. Aex ■ 63 " ©P "♦" ^P " ^s Case IV. AeX " Gg - Ag - Ap 



Case V. Aex - ^s - Gp 




LANETARY CAP 



Figure 12, Possible Penumbral Configurations 
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ONE llALF OF THE 
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ONE HALF OF THE 
SPHERICAL LUNE, Ap 

CIRCLE OF PROJECTED 
PLANETARY CAP 



CIRCLE OF 
PROJECTED 
SOLAR CAP 



GREAT CIRCLE 



Figure 13. Area of Spherical Lune, A 
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Aj), where Ag is the angular area of the smaller lune formed by the 
soiar cap and the great circle, and Ap is the angular area of the 
smaller lune formed by the planetary cap and the great circle. 

Figure 13 illustrates the geometry and the variables used to 
compute Ag. The method of computing Ap is identical except that 
all computations are done in the (x', y' , z') coordinate system 
instead of the (x, y» z) coordinate system. 

The area of the surface ABC of Figure 13 is given by: 



area of ABC - 



*0 %W 

f 
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a d4) 1 a sin G dG, -' ^q < 7 (79 



To get the angular area, Ag, one must double the area of ABC and 
divide the result by a^. Thus, 



As - 2 



^(*) 

d^ I sin 6 dG (80) 

g;(*) 



or 



where 



{cos [Oa((())] - cos [Gb((())] d^ (81) 



4) is the angular displacement in the x-y plane, 
positive counterclockwise, 

e^Cf^)) is the angular displacement of side BC (an 
arc of the great circle) from the z-axis, 
positive counterclockwise, 

6^(4)) is the angular displacement of side AC (an 
arc of the small circle of the solar cap) 
from the z-axis, positive counterclockwise. 

Before Equation (81) can be evaluated, expressions for cos [G^(4))] 

and cos [eij((())] must be found. The equation of the email circle of 
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the solar cap Is: 



x2 + y2 - H (2a - H) . 



Using the following transformation: 

X ■ a sin cos (J) 

y ■ a sin sin i> 

2 ■ a cos 

Equation (82) becomes: 

a^ (sin^e cos2(^ + sin^G sin^^) - H (2a - H) 

or a^ ( 1 - cos^G) - H (2a - H) 

For the small circle of the solar cap, 9 - O^t ^ind Equation (84) 
becomes 



cos 01 



a - H 



(82) 



(83) 



(84) 



(85) 



(86) 



The equation of the great circle is 

(x")2 + y2 . a2 

and using Equation (53) for x" we obtain: 

(x sin Oc + 2 cos Gq)^ + y^ « a^, (87) 

From Equation (75) we note that everywhere on the great circle 

z tan 0G * x, (88) 

Hence, after substitution of Equation (88) into Equation (87) we have 

(2 tan 0G sin 0G + 2 cos 0g)^ + y^ - a^. (80) 

Inserting the polar coordinate relations from Equation (83), yields after 
some reduction: 

cos^G (1 - cos^Gq sin24)) - cos^Gq cos^* (90) 

or, since G - 0a(4>) along the great circle 
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COS 0«(4>) 



cos 6q cos (J> 
J \ - cos^Gq sin2<> 



(91) 



Now we can evaluate Ag of Equation (61) using the value of cos 9^ 



and cos Gu, above. 



As- 2 



*0- 



cos Gq cos 4> 



r- 



a - H 



cos^Gg 8in^cJ> 



d4> 



Rearranging we have: 

f cos ^ dji o ( r-i H-i , 
As - 2 I , ' ' ^ - 2 (1-7) d^ 



Letting sin $ - x, c - sec 6q 



/sec Sfj - sin^cji o 

■ sec 

sin 4>Q 



As - 2 



dx 



/^^T7^ 



-2 ri-ii] 



(Po 



2 sin 



.^ ' sin 4»o 

L 



sec Gq 



-2 (l-f) ,c 



(92) 



^s - 2 sin~l [sin 4)o cos %q\ - 2 (^.q (l - j) 



From Equation (82) and Figure 13, it is evident that 



xp 



cos 4-0 - / I ^ ■ /H (2a - H) 



/xp2 + yp 



yp 



sin 4>o 



yp 



yxp2 + yp2 /H (2a - H) 



, yp > 



(93) 



(94) 



(95) 
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tlierefore. 



<^0 - sin 



-1 



yp 



/H (2a - 11) 



1 *0 IT 



(9fO 



Using Equations (93), (95), and (96) the value of Ag is obtained. 
Notice that the above development was done under the restriction* 
^ i, ^0 « T* ^^ which case the entire lune of Intersection, A^, is in 
the first and fourth quadrants. If the intersection of the two cans 
occurs in the second quadrant (as in Case I of Firrure 12), then <Pq 
is in the second quadrant and the entire lune of intersection lies in 
the second and third quadrants. Instead of recomputing this case, 
note that (see Figure 13) if the lune of intersection. As, lies in 
the second and third quadrants one could rotate the x-y plane about 
the z-axis by 180**, resulting in the problem tliat has .lust been ana- 
lyzed. Therefore, Equations (93), (94), and (96) give the correct 
value of the angular area As whether <pQ actually lies in the first 
or second quadrant. Notice, however, that if ^ JL *f*0 JLT ^^^^ 
xp ' and cos (t>Q > 0; but if y ^ <J> ^ tt then xp < and 
cos ^0 ^ 0* Thus, cos <t*Q can be used tc discriminate between the 
two cases just mentioned, 

A completely parallel development is used in order to compute 
Ap. The major difference is that all calculations are done in the 
(x'y'z') coordinate system. The resulting equations are: 



Tjl 

Ap « 2 sin"^ [sin ({>o' cos Gg'] - 24)o' [l - "TJ 

a 



(97) 



cos <^Q 



Xp COS e^^ - (a - 11) sin 0^ 



Al' (2a - II') 




sin <t>n' - ^P 




/H' (2a - H') 





(98) 



(90) 
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■?o' 



sm 



yp 

/li' (2a - H') 
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wiiere the primed quantities are related to the planet and Of;* is 
i^iven by Equation (78), 

Finally, note that many quantities are expressed in terms of <i , 
the radius of the great sphere, and H or H' the depth of the solar 
or planetary cap. By use of Equations (60) and (61) we have 



a 



27r 



(J 00) 



( 1 D J ^ 



a 27r 
Thus ail quantities may be expressed in terms of '.;<, and 6p which 

are already known. This is most easily accomplished by settinp a ■= 1 

(since the radius ot the sphere is arbitrary) and letting 11 » — 
6p -^ 

and H* - — . 

The compulation of ^^ux* ^^^ exposed angular area of the solar 
cap is shown for each of the five cases of Figure 12 in Table 11. Note 
that the discriminnts used to determine which case is to be evaluated 
are yp t cos (Pq^ and cos i q\ 

Finally the penumbral illumination factor is given by Equation 
(37): 

Aex 



(JOJ) 
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NUMLKICAL INTEGllATION 

The ccjuatioas of motion for a space vehicle are second-order 
diifereniial ccjuaLions which describe the accelerations arising; i ror 
the forces actin^^ on the vehicle. Accounting*, for the primary ijra\ i- 
rational field of the reference body and four types of perturbative 
acct lerations, the equations of motion (from Equation (2)), 

^^ - - '^ + ?1 + ?2 + P3 + ^ (2) 

If the perturbations are considered to be zero, the ri£;ht hand 
.T.ice of Equation (2) reduces to one term and the vehicle will follo^^ 
a Rcplerian orbit which may be described In closed forn In terras of 
it^ true or eccentric anomaly. Usually, however, part or all ol the 
perturbations are included and Equation (2) must he numerically inte- 
[;ratcd , 

There arc two basic methods by which the integration of Equation 
(2) may be formulated* Encke's .Tiethod and Cowell's method. If Equa- 
Lion (2) were to be numerically integrated in a straii^ht-f orward 
mjnncr» the integration would be knoi>m as Cowell's method. The sir- 
plicity of this method is offset by the larfe accelerations which must 
be integrated. As a consequence of the acceleration ma^.nitudes, small 
time increments have to be used in the Integration, and machine round- 
off error accumulates rapidly. Independent evaluations at many cor- 

panies and universities have shown that Cowell's method requires more 

u.ii. r>ni[W»PP- 228-242] 
machine time than other perturbational schemes. See Baker -^ 

tor a further discussion of Cowell's advantages and disadvantages. 

Despite its drawbacks, Cowell's method is still widely used and 

is especially suited if the accelerations due to perturbations are 

as large or larger than the term due to the central force field (e.f., 

during reentry). The program contains both methods and the option is 

left to the user. 

Historically, Encke's method is older than Cowell's although 

the former is more sophisticated. Cowell's method requires a modern 
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high-speed computer to be practical ^ whereas P'ruke's was dovcjopctl 
for hand computation. In Encke's method, it Is assumed that the per- 
turbative accelerations, P^, are small compared to the reference 
body acceleration. Consequently, when drap, accclerntions are small, 
rhe solution of Equation (1) is a good approximation to the true orbit. 
Under these conditions, it is only necessary to integrate the differ- 
ence between the accelerations on the two-body orbit and the total 
accelerations acting; on the vehicle. The cquation^v of motion then 
become second-order differential equations descrlhinr the acceleration 
differences. Let 

C - R - Rtb (-LOB) 

where Rjj^ is the position of the vehicle in termvS of the two-body 
orbit. Then, 



R ^Tb" 



1^ 

+ I h 

i-l 



(104) 



Equation (lOA) is integrated to obtain C and ^, These qunntititc 
are then added to Rtb and RtB» respectively, to obtain the instan- 
taneous position (R) and velocity (R) of the vehicle. The quan- 
tity i is commonly referred to as the "Encke" term, 

Encke's Method 



Rewriting Equation (lOA) we have. 



"^TB 



- -► •* /RTB^ 

Rtb + ^ - ^ - R ^ — ^ 



3 
B 



+ I Pi 
i 



R^3 



(i-^)^-^]-r^i 



The above equation could be integrated directly; however, the 

term ^1 - '"7^) is not well suited for the nurerical crl culatior. 

3 
since Rt^/R^ is very close to 1, Hence we rewrite this terr hv 

setting 
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r *• 



- i 



R • R + (R • R^g - R • Rj^) - Rjp^ • K'j.j^ 

"TB ^B 

R . (Rtb + I^) - Rtb (^Tr> + ^) 4 • (2Rtb + 7) 
. ^-5 . ^ ^ 



'TB 



Then, (1 + q) ^ •= RVR'IB and 



^TB [(1+q) ^ + 1] 1 -^ Or) ^2 

^xpar.tlinp, c'..l nunarator 

R^ , 3q -^ 3q^ ^ q' 

"1 — ■'■ * — vr~ • 

^B 1 -^ (l-^'^O 

In bumniary, LiioKe's meLhod computes the deviations from the nomin^.l 
orbit by integrating 



V 



^i 



where 



TSq + 3q2 + q3^ ^ 
L 1 + (]+n) 2 J 

Deterir.ination of Kenlerian Orbit Vectors 

Lncke's laethod requires the calculation of the position and velo- 
city vectors, Rtb and Rtb» respectively, of the norinal Keplerian 
orbit. The position and velocity at time t may be written in terrs 
of tne position and velocity at an earlier time Cq as follows: 



(H)^) 



(106) 



and 



RTB(t) "• f RTB(to) + g RTB(to) 
RTB(t) - i. RTB(to) + 8 RTB(to) 



(107) 
(108) 



where f and g are explicit functions of the differential eccentric 
anomaly of the Keplerian orbit. The program uses Herrick's method to 
determine f, {, g, and g, and then computes RtbCc) and RTB(t) 
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from Equations (107) and (108). The equations for computing f, f, g, 
and g are given below; their derivations may be found in Battin'" -' 
or Baker^lO]. 

f - 1 - C/ro (109) 

Tq - |RTB(to)l (110) 

M - U 

• - /y" s 

A rg (111) 

g - 1 - C/a (112) 



where 



where 



M - /T" (t-to) - ro X + do C + coU (113) 

(llA) 
(115) 
(116) 
(117) 
A - ro + doS + coC (118) 

CO . ^ - 1 (119) 

i - ^ - ^ (120) 

a ro y 



'=->'^ l/l - ',a * 6717 - 8^ -^ • . 


. .] 


"■>^''3' -5U-^7!a2-9fcT-^- ■ 


. .1 


a 

• 

RTB(to) ' ^B(to) 
do - " J=^ 





a is the semi-major axis of the Keplerian orbitt 
y is the universal gravitational constant, 
and 
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vo - l^iiCto)!. 

Solution of Equations (109) through (112) rsqulras th« datarul- 
nation of X. Equation (113) la aolvad for X by vrltlng 

F(X) - M - roX - dflC - cqU • 
and 

r'(X) - - ro - doC - cqU' • 0. 

Salactlng an Initial ••tlHat* of X to b« X|, tha root of 
F(X) - may be datarmlnad by a Nawton-Raphaon Itaratlvc procaaa, 
l<e., 

Xl^-i • Xi - P<^l>/7'(Xi) 

From Equations (114) and (115) 

C - X - U/a 

U' • C, 



Slnca Equations (114) and (115) ara aarlaa axpanalona Is powara 

7. a. 

an ellipsa, 



X2 
of -^, « method muse b« found to limit tha alza of thia tavii« For 



X - (E - Eq) /r 

vhera £, and Eq ara tha accantrlc anomaliaa at tlmaa t. ami to 
reapactivaly. Hcnca, by updating tha aj^^ch to at fraquant intarvals. 

tha diffarance £ - Eq ia kapt small* By similar raaaoniagt fra* 

X^ 

quant updating will kaap ^^* amall for hyparbolic orbita. 

Cowslips Method 

As described bafore« tha general equatlona of metion ef a epace 
vehicle are: ^ 

i - - ^4^ I ^i <121) 

^ i-1 
In Covall*8 method, theaa equatlaas are Intagratad* ualng nuaeri* 

cal techniquea* to obtain the Iii8taataa4i0us poaltlen and velocity of 

tha vehicle. 
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The program performs the integration using the same techniques as it 
does for the Encke method. The Runge-Kutta starting procedure pro- 
vides the initial data, and the Nordseick method is used as the long- 
term integration procedure. 

Integration Technique 

Equation (2) in Cowell's method or Equation (105) in Encke 's 
method must be numerically integrated by the program. This process 
is divided into two stages, a starting procedure and a long-term pro- 
cedure. The long'^term numerical integration procedure requires know- 
ledge of previous data points. Thus, the starting procedure is needed 
to provide the initial data points for the long-terra procedure. 

The procedure chosen by Sperry-Rand for long-term integration was 
based on their experience with methods used in the ITEM and MINIVAR 
programs, and the results of comparative testing. The justification 
of their choice is presented below as it appears in the original docu- 
ment on the SPACE program^ -^ . 

'*The long-term numerical integration procedure presently in use in 
the ITEM and MINIVAR programs is an Adams sixth-order predictor method 
(without corrector) for second-order differential equations. It was 
desired, however, to test a broader class of procedures before deciding 
on one for use as the long^^term numerical integration procedure to be 
used in the program. Accordingly, a program was written to test pre- 
dictor, predictor-corrector, and iterated predictor-corrector, i.e., 
repeated application of correctors, techniques of various orders of 
approximation and with or without modifiers. 
'*rhe following results were obtained: 

(a) Modifiers were found to leave the error unchanged. There 
is, therefore, no reason to use them. 

(b) As time interval increases, there is more tendency for the 
solution to become unstable. Error increases with a large 
power of the time interval. 

(c) As the degree of the approximating polynomial increases, a 
decrease in stability is noted. Error decreases about 3:1 
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for unit increase in degree of the approximating polynomial. 

(d) Predictor-only methods (degree and time interval held fixed) 
are about 30:1 less accurate than predictor-corrector meth- 
ods, i.e.p one application of the corrector removes 97% of 
the error in the predictor. A second iteration of the cor- 
rector does not reduce the error, but does Improve stabil- 
ity at the expense of a 2:1 increase in running time. The 
increase in running time is intolerable; therefore, a pre- 
dictor-corrector method will be used with no iteration of 
the corrector. 

(e) Either a fifth or sixth degree approximating polynomial 
yields a good compromise between accuracy and stability. 

**The final choice of a long-term integrating procedure involves 

considerations other than only accuracy and stability: 

(a) The ease of transforming the output of the starting proce- 
dure to the form of input starting data needed by the long- 
term procedure. 

(b) Whether the long-term procedure can easily accomodate a 
change in the time interval. 

(c) Whether the long-term procedure can easily interpolate 
to find conditions at an intermediate time at which data 
are desired. 

"There are at least three forms in which the Adams long-unrm pre- 
dictor-corrector formulas can be written. The only difference is in 
mathematical form; therefore, the accuracy and stability are the same 
for all three. 

"The first form is the conventional one in terms of the successive 
backward differences; the second is in terms of the successive values 
of the function. The third is due to Nordsieck and uses the succes- 
sive higher derivatives of the approximating polynomial. Each of the 
three forms has certain advantages and certain disadvantages which 
will be discussed now. 

"The backward difference form is fairly easy to start but inter- 
polation is somewhat difficult, and it is virtually impossible to 
change intervals except by use of the starting procedure. The suc- 
cessive value form of the method is trivial to start up, but inter- 
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polation involves the Lagrangian interpolation formulas, and changing 
intervals is, again, almost impossible. The inability to change in- 
tervals immediately after starting causes, as in the present ITEM pro- 
gram, a situation where the starting solution is called 28 times, but 
used only 7 times. 

"The Nordsieck method is fairly difficult to start, but very amen- 
able to arbitrary changes of time intervals and to interpolation to 
intermediate points. Five points are all that are needed to start 
after a change in time interval (of about 4:1). Due to its versatility, 
the Nordsieck method of degree 5 (called m " 6 by Nordsieck) without 
iteration and without choice of interval is used in the program." 

Starting Procedure (Runge-Kutta-Gill, RKG, Method) 

The program uses the Gill modification of Runge-Kutta*" -I 
for the starting procedure. Runge-Kutta methods are widely 
used and the differences between Gill's method and others are minor. 
Therefore, only the formulas as given by Gill and not the deriva- 
tions are presented. The method was chosen because it introduces some 
simplicity and error reduction over similar methods. 
Given the differential equation 

g-f(x. y) 
and the initial values xq and y(xo), y(xo+h) is calculated by 

y(xo+h) - y(xo) + ^ ko + ^ [1 - J\\ k, + ^ [1 + J\\ k2 + ^ k3 , (122) 
where 

ko - h f(xo, yo) yo ■ y(xo) 

ki - h f(xo + Y* yi) yi - yo + "2 ko 

k2 - h f (xo + \, y2) y2 - yo + [- 7 + f\\ ko + [1 - f\\ ki 

k3 - h f(xo + h, yg) ys - yo + [- f\^ ki + [1 + J\\ V^ (123) 
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The above procedure could be used directly to compute y(x(j + h) > 
y(xo + 2h)» ..» J etc. However^ Gill has introduced a "bridging" 
technique between one entry and the next. This modification increases 
the accuracy with little increase in complexity. The new procedure*" -^ 
may be summarized as : 

Given xo» yo» y' (xi) - f(x^, y(x^)) 

calculate for i « 0» 1» 2» 3 

Ki - a^(y' (xi) - bi qi) 

^i+1 - Xj_ + Ti h 

yi+1 - y(5Ci+l) - yi + h Ki 

qi+1 " qi + 3 Ki - ci y'(xi) 



where 

1 
2 



ao " T bo - 2 



ai - 1 - /^ bi - 1 



/I 



Co 


1 
■ 2 




To ' 


_1 


ci 


- 1 - 


■/T 


Tl 


- 


C2 


- 1 H 


^/T 


T2 


1 

2 


C3 


1 




T3 


- 



a2"l+/2' b2-l 

1 , 

a3 - ^ b3 

For the first entry into the procedure^ qg " 0> and the proce- 
dure is the same as given in Equations (122) and (123). For subse- 
quent entries » qg is set equal to the previous q^. 

It should be noted that the program must apply the method to the 
three second-order differential equations given by Equation (105) 
(Encke's method)^ or by Equation (121) (Cowell's method). The program 
treats the problem as one of six first order differential equations 
given by: 

yi - fi(t» yu y2» •••» ye) i - 1» 2 6 (124) 

where 
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yi(t) - x(t) y^(t) - x(t) 

y2(t) - y(t) y5(t) - y(t) 

y3(t) - z(t) y6(t) - i(t) 

where (x, y, z) and (x, y, z) are the geocentric position and velo- 
city vectors (in Cowell's method) or the perturbations from the nomi- 
nal orbit (in Encke's method). 

The procedure is programmed as follows: 

1. qi -^ i - 1, 2, ..., 6 

2. y - (yi, y2. •••» ye) - (R(to), R(to)) 

3. For k - 1, 2, 3, 4 

• • •• 

A. y - (yi. ya. •••• ye) (R. R) 

B. For i - 1, 2, ..., 6 

1. K -^ ak (yi - bk • qi) 

2. yi "^ yi + K . h 

3. q^ "^ qi + 3 • K - ck • yi 
c. R "^ (yi, y2. ya) 

• 

R *■ (yi*. ys. ye) 
t -^ t + Tk • h 

call DERIV, i.e., cal. R - F(R, R, t) 

• •• 

Hence the program exits with R(to + h) » R(to + h) , and R(to + h) . 

Also, steps 1 and 2 are only executed upon the first entry into the 

procedure. For second or later entries, the initial value of q^ is 

-► 
the value computed during the last entry for K ■ A. Clearly, y 

-► -► 
already contains (R(t), R(t)) and need not be reset. 

Transition from Starting Procedure to Long«*Term Integration 

The starting procedure yields the solutions of the six differen- 
tial equations and their rates of change at six successive times. It 
is necessary to transform these data into the form required by the 
Nordsieck long-term numerical integration procedure. 
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For each first order differential equation, the Nordtieck. method 
requires the following five higher derivatives evaluated at t - tg: 

h y(to) 



a(to) 
b(to) - 
c(to) - 
d(to) - 
e(to) - 



2! 
h2y(III) (to) 



3! 
h3 ydV) (to) 



A! 

h^ y(v) (to) 
5! 

h5y(VI) (to) 
6! 






J 



The RKG starting method provides data for y(t) and y(t) at 
the six time intervals up to and including to» i«e«» RKG provides: 

y(to) 



y(to) 
y(to - h) 
y(to - 2h) 
y(to - 3h) 
y(to - Ah) 
y(to - 5h) 



y(to - h) 
y(to - 2h) 
y(to - 3h) 
y(to - Ah) 
y(to - 5h) 



The required values for a(to), b(to), c(to)» d(to)» and e(to) 
will be found by using Lagrange's Interpolation formula to fit a 
power series of degree five to the y(t) data provided by the RKG 
method. The power series will then be successively differentiated to 
obtain the required data. Let 

X - t - to 
h 

and let primes denote derivatives with respect to x. Therefore, 



(125) 
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^'(t) - h y(t) 

y"(t) - h2 y(t) 

y"'(t) - h3 y(^^^(t) 

y""(t) - h** y(V)(t) 

y"'"(t) - h5 y(VI)(t) 

y"""(t) - h^ y^^^^^(t) 

From Lagrange's Interpolation Formula, 

5 
y(x) - I Fi(x) yi 

i-0 






where 



Fo(x) 



(x+1) (x+2) (x+3) (x+A) (x+5) 



120 



F, (x) - - x(x^2)(x+3)(x^-A)(x^-5) 



FjCx) 



F,(x) 



x(x+l)(x+3)(x+A)(x+5) 
12 

- x(x+l)(x+2)(x+A)(x+5) 
12 



V r^^ - x(x+l)(x^•2)(x^•3)(x+5) 



(126) 



(127) 



(128) 



2A 

F (x) - - x(x+l)(x^•2)(x+3)(x^•A) 
^ 120 

and yi are the values determined by the RKG method. Multiplying 
out the factors in Equation (128) yields 
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Fn(x) - x^ + 15 x** + 85 x^ + 225 x^ + 27A x + 120 

Fi (X) - - U' -^ 1^ P^** -^ 7^ ?^^ •*• ;?^ IS^ + 120^) 

F2(x) - x5 + 13 x** + 59 x^ + 107 x^ + 60x 

F3(x) . - ^?^' + 12^"* + ^1 x3 ■*■ 7? ?^^ ■*■ ^f}^) 

P^(^) . ^5 + ^ x-* + 4; x3 + 6 p^2 ^. ^0^ 

t, / X - (x^ + 10 x** + 35 x^ + 50 x^ + 24x) 

F5(x) ^ ^ ; 



(129) 



From Equations (125) and (126) 



a(to) 



iUloi 



2! 



b(t.) - ^ 



c(to) - 
d(to) - ^ 



v"' (0) 



(130) 



e(to) - 



41 




,1111 


CO}. 


5! 




JMM 


m 
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J 



-0- 

Successively differentiating Equation (127) with respect to x» 
(using Equation (129)), setting x - 0^ and substituting into Equa- 
tion (130) yields the following in matrix notation; 



2a(to) 




-24 


150 


-400 


600 


-600 


274 




y(to - 5h) 


3b(to) 




-50 


305 


-780 


1070 


-770 


225 




y(to - 4h) 


4c(to) 


. -L- . 


-35 


205 


-^90 


590 


-355 


85 




y(to - 3h) 
















t 




5d(to) 


120 


-10 


55 


-120 


130 


- 70 


15 




y(to - 2h)| 


6e(to) 




- 1 


5 


- 10 


10 


- .5 


1 




y(to - h) 
y(to) 



or 
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a(to) 
b(to) 
c(to) 
d(to) 
e(to) 





-144 


900 


-2400 


3600 


-3600 


1644 






-200 


1220 


-3120 


4280 


-3080 


900 




1 


-105 


615 


-1470 


1770 


-1065 


255 
















« 


1440 


- 24 


132 


- 288 


312 


- 168 


36 






- 2 


10 


- 20 


20 


- 10 


2 





















y(to - 5h) 

y(to - 4h) 

y(to - 3h) 

y(to - 2h) 

y(to - h) 

y(to) (131) 



Equation (131) is used to compute the sets of coefficients [a(to)» 

T 
b(to)f x(to)» d(to)» e(to)] for the six differential equations given 

by Equation (124). 

Long"-Tenn Integration (Nordsieck Method) 

The Nordsieck method is used to continue the solution of Equation (124) 



dyi 

"J^ " fi(t, yi, ...» ye) l " 1. 2, ..., 6 

once begun by the RKG method. 

Considering the single equation 



(132) 



dt 



f(t, y) 



and approximating its solution by a polynomial of degree five 
yP(to+h) - y(to) + h[y(to) + ^ y(to) + ~ y(l")(to) 
+ Try(^^>(t:o) +^y(v)(to) +^y(vi)(to)] 



A! 



5! 



6! 



(133) 
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where 

h is the integration step size, 
to is the value of t at the last integration* 
Substituting from Equation (125) 

y^(to+h) - y(to) + h[f(to» y) + a(to) + b(to) + c(to) + d(to) + e(to)] (134) 

Differentiating Equation (133) with respect to t = t^. + h and using 
Equation (125) 

f^ - f (to» y) + 2a(to) + 3b(to) + 4c(to) + 5e(to) (135) 

Having obtained a predicted value of f(to + h, y(to + h)) a 
value of f(to + h) ■ f (to + h, y^) is computed and combined with 
f^ to give the Nordsieck corrector, 

Ki h[f (to + h) - fP] 
where 

^1 -ii?- 0-315591931 

The values of a(to + h) , b(to + h) , c(to + h) , d(to + h) and 
e(to + h) are then computed in terms of their values at t - to 
by 

a(to + h) - a(to) + 3b(to) + 6c(to) + lOd(to) + 15e(to) 

+ K2[f(to + h) - fP] 
b(to + h) - b(to) + 4c(to) + lOd(to) + 20e (to) 

+ K3[f(to + h) - fP] 
c(to + h) - c(to) + 5d(to) + 15e(to) + Kjf(to + h) - fP] 
d(to + h) - d(to) + 6e(to) + K5[f(to + h) - fP] 
e(to + h) - e(to) + K6[f(to + h) - fP] (136) 



where 












K2 


- 


137 
120 


- 


1. 


.141666667 


K3 


- 


5 
8 


- 


0. 


.625 


K^ 


- 


17, 
96 


- 


0. 


.1770833333 
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Ks --^- 0.025 

Ke - rjj^ - 0.00138888889 

The Nordsleck method may now be summarized as follows. For each 
entry, y^ is computed by Equation (134), f by Equation (135), and 
f (to + h) " f (to + h, y^) by the functional relationship implied by 
Equation (132). Then y(tQ + h) is found by 

y(to + h) - yP + Ki h(f(to + h) - f P) . 

Finally the coefficients a(t) through e(t) are updated by Equa- 
tion (136). 

The integration interval, h, is readily changed; the change 
being accomplished by using new values of a(t), b(t), c(t), d(t), 
and e(t). These new values are obtained from the following equations: 

^ ho 
an(t) - Bao(t) 
bn(t) - B2bo(t) 
Cn(t) - B3co(t) 
dn(t) - B-doCt) 
en(t) - B^eoCt) 

where the subscripts n and stand for new and old, respectively. 
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COORDINATE SYSTEMS AND TRANSFORMATIONS 

The equations of motion of a vehicle are given by Equation (2)t 
It is necessary to express the vectors used in this equation in an 
inertial coordinate frame, i»e., a coordinate system in which inertial 
forces due to the system's acceleration are absent or at least negli- 
gible. A coordinate system rigidly attached to the Earth is inappro- 
priate, since such a system experiences noticeable accelerations due 
to the daily rotation of the Earth. 

The inertial frame chosen for this program is a so-called "mean 
equinox of base date" system which is determined by the vernal equi- 
nox of 0.0 hours 1 January of the year subsequent to the. epoch time 
(initial input time) used in a given run. This particular base date 
system has been chosen as a basis for calculation because the plane- 
tary, lunar, and solar coordinates are written on tapes in that coor- 
dinate system. Rather than transform the tape information, the vehi- 
cle initial conditions and the oblateness accelerations are transformed 
into the base date system. 

The definition of the base date system and various auxxliary 
Earth referenced systems is given below along with the transformations 
between them. 



Terminolo 



sx. 



First, consider the celestial sphere (Figure 14) which is a 
sphere of infinite radius with the Earth at its center upon which the 
positions of stars, planets, the Sun, and the Moon are projected. It 
is sometimes convenient to pick the center of the sphere to be at an 
observer or the center of the Sun, but for the purpose of our discus- 
sion its center will always be taken to be at the geocenter. The 
stars may be taken to be fixed on the sphere (their small proper 
motions being neglected) while the Sun, planets, and the Moon appear 
to move and describe certain paths along its inner surface. The ap- 
parent path of the Sun during the course of a year on the celestial 
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sphere is a great circle called the ecliptic . The north celestial 
pole is the point where the extension of the Earth's north pole inter- 
sects the celestial sphere and the south celestial pole is defined 
similarly. The great circle where the extension of the Earth's equa- 
tor intersects the celestial sphere is called the celestial equator . 
The point on the celestial equator where the apparent path of the sun 
crosses into the northern half of the celestial sphere at the begin- 
ning of spring is called the vernal equinox . This is one of two 
points where the ecliptic and equator intersect. In the following 
discussion the word celestial will usually be dropped and we will 
just speak of the north pole or the equator. 

The equator and ecliptic are constantly in motion due to the ef- 
fects of nutation and precession. It should be pointed out that this 
is astronomical precession which is distinct from the force free pre- 
cession due to the flattening at the poles. The former is caused by 
the net torque on the equatorial "bulges" due to the gravitational 
attraction of the sun and moon. The torque is quite small so that the 
precession is extremely slow - the period being 26,000 years compared 
to the rotational period of one day. The total applied torque is not 
constant in time, because the torques of the Sun and Moon have slightly 
different directions to the ecliptic and vary as the three bodies move 
around each other. As a result, there are irregularities in the pre- 
cession, commonly designated as astronomical nutation . The astronomi- 
cal nutation is not to be confused with the "true nutation" which is 
present even in the force-free precession of the Earth's rotation 
axis. In the subsequent discussion the astronomical precession and 
nutation will simply be referred to as precession and nutation. 

Precession may be separated into luni-solar precession, and 
planetary precession. The luni-solar precession causes the north 
pole to move around the ecliptic pole . The ecliptic pole is the 
point of intersection of the extension of the normal to the ecliptic 
plane at the geocenter, and the celestial sphere. Nutation is an 
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Figure 14. Celestial Sphere 
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Figure 15. Precession 



irregular motion of the pole with a period of about 18.6 years. These 
two effects combine to model the motion of the true equator . The ac- 
tual equator at any time is callad the true equator of date. The word 
date here refers to any arbitrary time. The mean equator of d«C« i« 
defined by the location of the equator when the effects of nutation 
are removed (i.e.| where the true equator of date would be if there 
were no nutation)* 

The gravitational action of the planets on the Earth causes the 
plane of the Earth* s orbit to slowly precess. This motion appears 
from the earth as a slow precession of the ecliptic, and is called 
planetary precession. This effect causes the obliquity (the angle 
between the equator and the ecliptic) to decrease at the rate of 
about 47" a century. Therefore, to be precise we refer to the eclip- 
tic of a specified date. 

The effect of precession is illustrated in Figure 15. The mean 
pole, mean equator, and ecliptic at time to are labeled base date 
and at a later time, t, are labeled date. Luni-solar precession 
causes the mean equator to move between base date and date, while 
planetary precession moves the ecliptic. Figure 16 illustrates how 
nutation causes the true equator of date to differ from the mean 
equator of date. 

In this discussion we will use several coordinate systems which 
will be defined in terms of the t ^ifue equator and ecliptic and mean 
equator and ecliptic at a given time. The true vernal equinox (or 
true equinox) is defined by the point of intersection of the true 
equator and ecliptic. The mean vernal equinox (or mean equinox) is 
defined by the point of intersection of the mean equator and the 
ecliptic. The true equinox and the mean equinox are points on the 
celestial sphere that experience a slow motion and hence we must 
specify a time to indicate their exact positions. 
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Coordinate Systems 

The five coordinate systems used by the program for generating 

trajectories are: 

^^^ '^^^ ^^^^ ^<l^^^ox of base date system . This coordinate 
system Is tVie inertia! system emplcyed by the program 
and employs unit basis vectors x, y, and z defined 
as follows: 

X is a unit vector directed towards the mean vernal 
equinox of base date, 

z is a unit vector normal to the mean equatorial 
plane of base date, positive in the northern 

hemisphere. 

y is a unit vector orthogonal to x and z, 

^^^ The mean equinox of date system. This system employs 

unit basis vectors xj,^, y>t» and zm defined as 
follows: 

X}] is a unit vector directed toward the mean 
vernal equinox of date, 

zil is a unit vector normal to the mean equatorial 
plane of date, 

yj.| is a unit vector orthogonal to x^ and z]kf. 

^^^ "^^^ true equinox of date system . This system 

employs unit basis vectors x-p, yx» and zx 
defined as follows: 

XX is a unit vector directed toward the true 
vernal equinox of date, 

Zj is a unit vector normal to the true equatorial 
plane of date, 

y-p is a unit vector orthogonal to xj and zx« 

(4) The Greenwich coordinate system . This system 

employs unit basis vectors xG» yCt and zq 
defined as follows: 

XQ is a unit vector directed toward the inter- 
section of the Greenwich meridian with the 
true equator of date, 

ZQ is a unit vector normal to the true equator 
of date. 
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yc is a unit vector orthogonal to xq and zq. 

(5) The topocentrlc coordinate system . This system employs 

three orthogonal basis vectors Xg^ y^ and z„^ 
directed toward the east^ norths and local vertical^ 
respectively > at a point on the Earth's surface 
(considered as an ellipsoid), it Is only used for 
observation calculations and Its details are dis- 
cussed under OBSERVATIONS. 

The differences between the above coordinate systems are due to 
dynamical effects of the Earth's motion or to the geometry of the 
Earth,' It Is possible to express a 3 >^ 1 vector of position or velo- 
city written In the basls^ of one coordinate system In terms of another 
coordinate basis by applying 3 >^ 3 orthogonal matrix transformations. 
Table III below lists the transformations used by the program. R is a 
3 >< 1 vector whose subscript indicates the basis in which the vector 
is expressed* 

Table III. 
Transformation Matrices. 
MATRIX EFFECT 



[G] Earth Geometry 

[y] 



FROM 


TO 


EQUATION 


h ys ^s 


xG yc ZG 


Rg - [G] Rg 


XG yc ZG 


XX yi zx 


^x - [y] Rg 


^T yx ^x 

XM yw ZM 


XM yM ZM 
X y z 


^M - [N] Rt 
R - [P] Rm 



Right Ascension of 
Greenwich from True 
Equinox of Date 
(Daily Rotation) 

[N] Nutation 

[P] Precession 

The derivations of [y], [N] ^ and [P] are given below. It 
should be recognized that these three matrices are functions of time. 

Sidereal Time 



When computing the inertial coordinates of a point or station 
on the surface of the rotating Earth at a particular time> a relation- 
ship must be found between the Greenwich meridian (the great circle 
that passes through the true north pole and Greenwich) and the vernal 
equinox at that instant. The angle measured in the equatorial plane 
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from the Greenwich meridian to the true vernal equinox is called the 
apparent Greenwich sidereal time . 

The right ascension of the mean equinox measured from the true 
equinox is called the equation of the equinoxes or nutation in right 
ascension and is usually denoted by 6a, The relationship between 
the apparent sidereal time and the mean sidereal time is: 

Mean sidereal time - Apparent sidereal time - 6a, 
6a is always measured along the true equator from the true equinox 

eastward. The apparent and mean Greenwich sidereal times are tabu- 

[41 u 
lated in the American Ephemeris and Nautical Almanac -* for O'^ of 

each day of the year along with the equation of the equinoxes, 6a. 

Each entry is given in hours, minutes and seconds of time to 0^.001, 

where 

1^ - 15\ 1^ - 15', and 1^ - 15", 

Apparent Greenwich sidereal time, Oq, since it is measured 
from the true equinox, is affected by nutation and hence changes at 
an irregular rate (see Figure 17), Since mean sidereal time, Sg, 
is measured from the mean equinox, which is affected only by preces- 
sion, it increases at almost a constant rate with only a small accel- 
eration and may be calculated at 0^ of any day by the formula 

Sg - 6^ 38m 45!836 + 8,640,184!542T + 0?0929t2 (137) 

where T is the number of Julian centuries (36,525 days) from noon 
January, 1900 (Julian day 2,415,020,), 
An equivalent formula in degrees is 

eG(d) - 100:0755426 + 0!985647346 d + 2^9015 x 10"! 3 d^ (138) 

where d is the integer number of days past 0^ 1 January, 1950, 
Equations (137) and (138) may be reformulated as a function of the 
time past any conveniently chosen epoch resulting in an equation of 
the same form with different coefficients. 

Evaluating the above formul* for eG(d + 1) and 6G(d), we 
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can find the angle through which the Greenwich meridian rotates in a 
day with respect to the moving mean equinox. 

eG(d + 1) - 0G(d) + 360: - 360:985647346 + 5? x lO"^^ ^^ (139) 

in which the additive term 2?9015 x 10"^^ has been neglected. Di- 
viding Equation (139) by the number of seconds in a day (86,400) and 
neglecting the final term whose contribution is insignificant over a 
period of only a few years ( < 25 yrs,), we have Wg, the rotational 
rate of the Earth with respect to the moving mean equinox, 

ws - 0;0041780 74622 sec."^ (140) 

Therefore, to find the mean Greenwich sidereal time at time 
d + T we use 

eG(d + t) - eG(d) + Ws T (141) 

where d represents the integer number of days past 0^ 1 January, 
1950, and t the number of seconds that have elapsed since of 
the day. 

The apparent Greenwich sidereal time 6q can be found at time 
t by 

^qM - eG(t) + 6a. (142) 

Now, big will be referred to as the sidereal rotational rate 
of the Earth. The word sidereal is used since the rate is with res- 
pect to the mean equinox, which is the reference point used for 
sidereal time. 

Transformations 

The effects of nutation and precession are fairly small and 
over a period of a few days the errors introduced by neglecting these 
effects is on the order of ten meters at most. Therefore, for short 
trajectory determination runs the program has the option of not in- 
cluding them by setting the [P] and [N] matrices equal to the 
identity matrix, [I]. The [y] matrix which includes the effects 
of daily rotation must be included, however, 
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. Figure 18. Formation of Nutation Matrix 
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Referring to Figure 17 the transformation from the Greenwich co- 
ordinate system to the true equinox of date system consists of a ro- 
tation along the true equator through the angle S^Ct) where 0G(t) 
is given by Equation (142). Thus the [y] matrix is given by: 



[y] 



If the effects of nutation and precession are to be neglected then 
egCt) replaces ©g^^^ ^^ Equation (143), and the reference system 
of the computation may be considered as the system defirted by the 
mean equinox of date (i«e«, the effects of nutation and precession 
are neglected) , 

Next, an expression for the nutation matrix [N] is found. 

In Figure 18, y is the true equinox of date and y is the 
mean equinox of date. The right ascension of the mean vernal equinox 
referred to the true equinox is 6a, the equation of the equinoxes. 
Notice that 6a is less than zero as it is drawn in Figure 18, since 
right ascension is measured positively toward the east. The nutation 
in longitude, Aij^, is the longitude of the mean equinox measured 
from the true equinox. Celestial longitudes are measured along the 
ecliptic, and since the positive direction is east, L^ is also neg- 
ative in Figure 18, The mean obliquity, e, is the inclination of 
the ecliptic to the mean equator, and e, the true obliquity, Is the 
inclination of the ecliptic to the true equator. The nutation in 
obliquity is 6e, where 

6e ■ e - e. 

From inspection of Figure 18 ^ and by using Napier's rules for 
right spherical triangles, 

sin(90** - e) - tan 6a tanCSO"* - L^) 
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(U3) 



or 



tan 6a ■ tan Ai)/ cos e 
Using the series expansion for the tangent function we have 
A ^ (<Sa)^ _^ 2 (6a) 5. (Ai^)3 ^ 2 (Ail;)^ ^ ^ 

Nowj 1 6a I, and |6i^| are less than li ^ lO""** radians. There- 
fore » the approximation 

6a - A\|; cos c 

should never produce an error greater than 1 ^ 10"^^ radians or 
2 X 10~^ seconds of arc. 

The nutation matrix [N]» is now derived which transforms co- 
ordinates referred to the true equinox of date into coordinates re- 
ferred to the mean equinox of date. This transformation consists of 
three rotation matrices. Using Figure 18, we note that the first, 
[A], is a rotation about the x-axis through the angle e. The sec- 
ond, [B], is a rotation about the z^-axis through AiJ;, and the 
third, [C], is a rotation about the x"-axis through - e. 
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- sin c 


COS c 
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[C] - 



10 
COS G - sin e 
sin c cos c 



[N] - [C][B][A] 



cos Ai^ 



sin h^ cos c 



sin Ai^ sin c 



1 


^4/ cos c 


Ai(/ sin e 


- L^ cos e 


1 


6c 


- Ai(/ sin e 


- 6e 


1 , 



- cos c sin h^ cos c cos Ail^ cos c cos c cos Ai{^ sin c 
+ sin G sin G - sin t cos t 

- sin G sin Ail^ sin g cos Atp cos g sin g cos Ail^ sin g 
- cos G sin t cos g cos t 

Since \&^l>\ < 10"** and |g - g| < 10""** the nutation matrix is 
often approximated by the following matrix: 



[N-] 



The error in the elements in [N'] should be less than one unit 
in the eighth significant figure'-^ » P- 3]^ 

Reference [l4] gives formulas for Ail^ and 6g which depend on 
the mean longitude of the Sun, the mean longitude of the perigee of 
the Sun, the mean longitude of the Moon, the mean longitude of the 
ascending node of the Moon, and the mean longitude of the lunar peri- 
gee. There are 69 terms in the expression for Ail^ and 40 terms in 
6g. Truncated expressions for A^ and 6e are given in Reference [l]. 

Finally, an expression for [Pj^ the precession matrix is de- 
rived below. Figure 19 Illustrates the effect of precession over the 
period of time t - tQ. In this figure, 90** - Co ^s the right ascen- 
sion of the ascending node of the mean equator at time t on the 
mean equator at to measured from the mean equinox of to; 90** + z 
is the right ascension of the node measured from the mean equinox of 
t; and 6 is the inclination of the mean equator at time t to the 
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Figure 19. Formation of Precession Matrix 
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mean equator of tQ. 

The precession matrix, [P], is now derived which transforms 
coordinates referred to yCt) to coordinates referred to Y(to). 

[P] is a result of three rotation matrices. From Figure 19 we 
see that the first, [A'], is a rotation of (90° + z) about the z- 
axls. The second, [B'], is a rotation about x through -9 and 
the third, [€']• is " rotation about z through Uq - 90°). 

[P] - [A'] [B'l [C], where 



[A'] 



cos (90° + z) sin (90° + z) 

sin (90° + z) cos (90° + z) 

1 



[B'l 



1 






cos (-9) 
- sin (-8) 




sin (-8) 
cos (-9) 



- sin z 

- cos z 




cos 8 
sin e 



cos z 

- sin z 

1 




- sin 9 
cos 9 



r 



[C] - 



COS Uo - 90°) 
- sin (Co - 90°) 




sin (Co - 90°) 
cos Uo - 90°) 




sin Co 
cos Co 




- cos Co 

sin Co 

1 



Hence, 



[P] - 



"- sin z sin ^o 

+ cos z cos 6 cos ^0 



cos z sin Co 

+ sin z cos 9 cos ^q 



- sin z cos Co c^s 2 cos Co 

- cos z cos 6 sin Co " ^^^ ^ ^^^ ^ ^^^ ^0 



- cos z sin 9 



- sin z sin 



sin 9 cos Co 
- sin 9 sin Co 
cos 6 



Reference [ 14] gives the following numerical expression for Co» ^» 
and 9 for the case when coordinates are precessed to a later epoch, 
Initial epoch, to - 1900,0 + To 
Final epoch, t - 1900,0 + to + t 

Co * (2304. "250 + l/'396 To) t + 0/'302 t^ + 0."018 t^ 
z - Co + 0."791 T^ 
9- (2004. "682 - 0."853 To) t - 0."426 t^ - 0."042 t^ 
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(145) 



where tq and t are in units of tropical centuries. 

When the initial epoch is 1900.0 + tq + x and the final epoch is 

1900*0 + Tq (i>e«y we are precessing a vector to an earlier epoch) 

one whould replace 
Co by - 2 
z by - Co 

e by - e 

in the above matrices [A' ] . [B'lt [C], and [P]» 

Combining the results of Equations (143), (144), and (145) with 
the relations in Table III, all necessary coordinate transformations may 
be performed by the program, 

OBSERVATIONS 

The SPACE-A program has a provision for a total of 25 observation 
types of which a number of ground-based observations are specified and 
described here. The specified observations include: 

azimuth, A 

elevation, E 

round-trip range, 2p 

one-way range rate, p 

hour angle, H 

declination, D 

il-direction cosine, i 

m-direction cosine, m 

X-angle, X 

Y-angle, Y 

range equivalent time. At 

range-rate equivalent time, At' 
The position and velocity of a vehicle is available at all times 
from the trajectory generation portion of the program. The station 
position and velocity is cpmputed from input data concerning its geo- 
detic coordinates. Therefore, one may consider the relative geometry 
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of the vehicle and station to obtain the observations. 

The program has the option to include the effects of refraction 
upon ground-based observations. A simple model of the Earth* s atmos* 
pheric effects upon observations is employed in which the index of re- 
fraction varies only with altitude. Neglecting the effects of refrac* 
tion causes the predicted elevation angle to be less than the observed 
elevation. Errors are also made in the computation of range and range- 
rate related observations. 

The method used to include the effects of refraction is to com- 
pute range and range rate from the relative geometry of the station 
and the actual vehicle position and velocity and then to. apply correc- 
tions based on the refraction model in order to obtain the range and 
range rate related observations. In the case of angle related mea- 
surements a slightly different approach is taken. A vector pointing 
from the station to the computed vehicle position is constructed. 
Then the refraction model is used to obtain a new vector pointing from 
the station to the observed vehicle position. All angle observations 
are then computed directly or indirectly using this new vector. 

Coordinates Used in Observation Calculations 

The input coordinates of the station (the station is assumed to 
be on an ellipsoidal Earth) are initially given as geodetic longitude, 
(Xq), geodetic latitude {^q) $ and height. Given this information 
the position and velocity of the station in the true system of date 
coordinate system is given by: 

Xj - [c + h] cos 4*0 cos (Xg + y) 

yx - [c + h] cos <\>G sin (Xq + y) 

zj - [cd-e^) + h] sin (^q 

XT - - Wg yx 

yx - uJe xx 

2X ■ 
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(1A6) 



c " ^e 



J \ - e^ sin^ 4>G (1A7) 

e2 - 2f - f2 (148) 



where 



^19 yT» ^T ^^c the true system of date coordinates of 
the station, 

^Gi ^G> ^ ^^^ t^^ geodetic coordinates of tha station, 

Y is the right ascension of the Greenwich 
meridian In the true system of date^ 

Re Is the equatorial radius of the Earth, 
and f is the flattening constant of the Earth- [ oqq on ) 

Using the precession matrix, [P] and the nutation matrix, [N], we 
can express the station coordinates In the inertlal (base date) coor- 
dinate system. From the position and velocity vectors of the station 
and the vehicle one can obtain the vectors of position and velocity of 
the vehicle with respect to the station expressed in the inertial sys- 
tem. Thus 

P - Pg - R (U9) 

• • • 

P - Ps - R (150) 

where p, p are vectors of the vehicle with respect to the 
station, 

->• ->• 

Pg, Pg are the vectors of the station, and 

R, R are the vectors of the vehicle. 
Note that all of these vectors are expressed in the inertial coordi- 
nate system. 

For observation calculations it is convenient to define a station 
oriented coordinate system which has a basis consisting of three unit 
orthogonal vectors pointed toward the east, north and local vertical. 
Using station oriented coordinates the unit vectors are: 



90 



xg - [1 0] 
^-[01 0]' 



2S 



-[0 1]\ 



east vector 
north vector 
up vector 



Expressed in terms of the inertial base date system these vectors arc: 



e 



[T] X, 



and 



east vector 
[T] yg north vector 
" - [T] zg up vector 
[T] - [P][N][Y][G][a] 



n 



(151) 
(152) 
(153) 
(154) 



where all vectors are 3x1 vectors. 

[a] is a 3 X 3 rotation matrix to account for 
misalignments between the station oriented 
coordinates and the true topocentric east, 
north, and vertical coordinates, 

[G] is a 3 X 3 matrix transforming coordinates 
expressed in the topocentric system into 
coordinates of the Greenwich system, 

[P],[N], and [y] are the 3^3 matrices of precession, nuta- 
tion, and right ascension of the Greenwich 
meridian which transform the coordinates into 
the inertial base date system* 



[G] - 



- sin Xq - sin ^q cos \q cos ^q cos Xq 
cos Xq - sin <(>G sin Xq cos (J)G sin Xq 
cos 4>G sii^ 4>G 



(155) 



where Xq and <(>g are, respectively, geodetic longitude and geo- 
detic latitude of the station. 

The final result of the above transformations is that the station 
oriented basis vectors (e, n, and u) are expressed in the inertial 
base date coordinate system. 

Computation of Observation Types 

In the following discussion all computations are made using vec- 
tors expressed in the inertial (base date) coordinate system. 

Azimuth and elevation measurements are illustrated in Figure 20. 
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Azimuth, A, is meaBured easterly from station north from ^ A ^ 2it< 
Elevation, £, is measured from the stations horizontal plant upward 
with a range " T i E i T« If P is the station to vehicle vector 
the azimuth and elevation are given by: 



A - tan 



Lp • nJ 



E " tan 



-1 



/ (P • e 



p * u 



)2 + (p • n)' 



(156) 



(157) 



If refraction is to be included in the computations, the vecter, p, 
is replaced by the vector p\ where p' is the vector from the 
station to the observed vehicle position as described below. 



u Local vertical 



y^^^ Vehicle 




^ n North 



e East 



Figure 20. Azimuth and Elevation 

^^® ^ound"trip range, 2e, is twice the distance fro« the station 
to the vehicle and its value is given byt 

2Q - 2|p| + 2Ap 

If the effects of refraction are not Included in the calculation 
2Ap - 0. If the refraction effects are Included Lp is a correction 



(158) 
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terra the value of which is given by Equation (202). 

The one-way range-yate ^ p^ is equal to the magnitude of the 
time derivative of the vector from the station to the vehicle and its 
value is given by; • 

P - ^TTr^+ ^P (159) 

If the effects of refraction are not included Ap - 0. If the effects 
of refraction are included Ap is a correction term the value of 
which is given by Equation (206) # 

The hour angle, H, is the angle between the station meridian 
and the projection on the true equator of the station to vehicle vec- 
tor measured in the true equatorial plane. It is measured westward 
from - IT ^ H ^ ir. The declination , D, is the angle made with the 
true equatorial plane by the station to vehicle vector. Declination 
is measured positive in the northern hemisphere and has a range 
" T — ^ — V Wour angle and declination are illustrated in Figure 21 
which shows the station at the center of a celestial sphere and the 
projected true equator on the celestial sphere. The hour angle and 
declination may be derived using spherical trigonometric relations in 
terms of geodetic latitude, azimuth and elevation (({) , A, E) , The 
values are given by: 



H - tan"^ 



sin A 



(160) 



cos A sin ({JQ - cos (jsq tan E 

D - sin"^ [sin <^q sin E + cos <^q cos E cos A] (161) 

The values of A and E used above are given by Equations (156) and 

(157). 

The JZ>-direction cosine , £, and m-direction cosine , m, are 

shown in Figure 22. The il-direction cosine is the cosine of the angle 

between tne station to vehicle vector p, and the station east vector, 

->- 

e. The m-direction cosine is the cosine of the angle between the sta- 

tion to vehicle vector p, and the station's north vector, n. The 
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values of l 



and m are given by: 

I - 



If refraction effects are included p' replaces p in the above 
equations f 

^^^ ^"^^1^ and Y-angle measurements are illustrated in Figure 23, 
The Y-dngle is the angle between the station to vehicle vector^ p^ 
and the perpendicular projection of this vector on the station's east- 
vertical plane. It is positively measured toward the norths negatively 
toward the south with limits - T — ^ ^1 T* "^^^ X-angle is measured 
between the vertical vector, u, and the perpendicular projection of 
the station to vehicle vector onto the station's east-vertical plane. 
It is measured from the positive vertical eastward with limit 
- IT j^ X ^ TT, The value of X and Y are given by: 



(162) 
(163) 




Y - tan 



If refraction effects are included p' replaces p in the above 
equations. 

The range time equivalent ^ At, and the range-rate time equiv- 
alents At', are included since the raw data from typical tracking 
systems are the time between a transmitted and received signal for 
range, and the time to count a given number of Doppler cycles for 
range-rate. In most systems, these quantities are first converted to 
range and range-rate units • However, it may be found useful in some 
cases to use the raw data equivalents* The values of At and At' 
are given by: 



(164) 



(165) 



(166) 
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Figure 21. Geometry Illustrating Hour Angle and Declination 
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Figure 22. i-cosine and m-cosine Figure 23. X-Angle and Y-Angle 
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Lt' 



1^21 - Ipll (167) 



where 



c is the velocity of light, 

jpil and IP2I are the ranges at the beginning and end 
of the measurements I respectively. 

Specific provision is made in the program for computing dopplcr 

shift frequency measurements. The equations for these observations 

are not included, since the method of computation can change depending 

on the specific hardware in use for a particular mission. 

Effects of Refraction 

Because of variations in the refractive index of the atmosphere 
the propagation path of a tracking beam is subject to refractive 
bending. 

Two models for the variations of refractive index are included 
in the program: one is a model for the troposphere, the other is a 
model for ionosphere. A numerical method due to Weisbrod and 
Anderson'- ■■ is employed in which the effects of refractive bending 
are determined by numerically integrating over the total propagation 
path, the index of refraction at each point being determined by the 
assumed model. The vector from the station to the current calculated 
position of the vehicle, p, is available from the program. Using 
numerical methods a correction to the elevation angle is found and a 
new vector, p', from the station toward the observed position of 
the vehicle is determined. 

In addition to angle corrections due to refractive bending, the 
effects on range measurements due to signal retardation, and the ef- 
fects on range-rate measurements due to refractive bending are also 
found • 

Index of Refraction Models 

In order to simplify computational problems, atmospheric models 
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representative of average conditions are employed in which the fol- 
lowing assumptions art made: 

(1) The gradient of the index of refraction varies only 
with altitude. 

(2) The index of refraction profile is approximated by a 
number of linear segments in which the length of each 
segment is very small compared to the Earth's radius, 

(3) The troposphere extends to 40 kilometers. 

(4) The region between the end of the troposphere and the 
beginning of the ionosphere is assumed to have zero 
refractivity (and therefore no bending or signal re- 
tardation occurs). 

(5) The ionosphere lies between a height ho (input data) 
and 2,000 kilometers. 

(6) The index of refraction is zero beyond 2,000 kilometers. 
In the tropospheric model, refractivity (N) is assumed to decay 

exponentially, with the ground index of refraction and the scale 
height as parameters. The equation for the tropospheric model is as 
follows: 

N - No e"^^" - (n - 1)10^ (168) 

where 



No is the refractivity at sea level, an input quantity 
for each station (usually 313), 

h is the height above the Earth, 

H is the scale height, an input quantity for each 
station (usually 7 kilometers), 

n is the index of refraction. 

For the tropospheric model, the refractive errors are considered 
to be independent of signal frequency since the index of refraction 
is virtually independent of frequency up to 30,000 MHz, 

In the ionospheric model, the index of refraction is dependent 
upon more parameters than those considered for the tropospheric model. 
The ionosphere consists of several belts of charged particles. The 
F layer is very much larger than any other layer, and contains a 
greater number of charged particles than the other layers. The F 
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layer is the one closest to the Earth's surface. It is subdivided 
into the Tl and F2 layers. In the ionospheric models the index 
cf refraction is primarily dependent upon the height, ho, of the 
base, of the ionosphere's F2 layer » the maximum electron density of 
the F2 layer, and the height of the maximum electron density of the 
F2 layer. 

Both index of refraction and the height, ho, are dependent upon 
diurnal, solar activity, seasonal, and geographical variations as well 
as other miscellaneous sporadic v^iriations. Unlike the tropospheric 
model, the refractive errors in the ionospheric model are frequency 
dependent. 

In constructing the model, the range of the signal frequency has 
been limited to frequencies above 100 MHz, since this range of spec- 
trxim both represents the situation of greatest interest and enables 
equation simplification. 

The relationship between the index of refraction (n) the angular 

frequency of the incident signal (w) , and the electron density in the 

ionosphere -' is given by j, 

1 '2 



n 






(169) 



Co ™ w^ 
where 

pg is the electron density per cubic meter, 

e is the electron charge (1.6 x 10"^^ Couloumbs), 

m is the electron mass (9.08 x 10"^^ kilograms) 

eo is the permittivity of free space (8.854 ^ 10"^^ 
farads/meter) . 

Using the first two terms of the binomial expansion as an approx- 
imation, the equation for the index of refraction reduces to: 

n - 1 - 40.3 ^ 10-12 (170) 

where f - ^ . Note that f has the dimensions of MHz. This 
equation is true for frequencies above the critical frequency, f^,. 
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which is defined as: 

fc - 8,97 Po ^ « 10"^ (MHz) (171) 

where pq is the BAximum electron density per cubic meter. 

From the definition of N in Equation (168), Equation (170) can 
be written as 

N - - 4,03 -| X 10*5 (172) 

The modal selected for electron density versus height consists 
of a parabolic variation below the height of maximum electron density 
matched to a hyperbolic secant profile above the maximum. The rela- 
tionships are as follows t 

Pe ■ Po [1 " (1 " ^)^1 ^0^ — ^ — 

p^ ■ sech T (<^ " 1) for o ^ 1 

where 

h - ho 

- 



h is the height above the Earthy* 

ho is the height of the base of the F2 layer (input 
quantity) » 

h|Q is the height of the maximum electron density in 
the F2 layer (input quantity). 

Figure 24 is a plot of the ionosphere model normalized with res- 
pect to o and 1/2 (Pg/Po)* "^^^ ^» ^n» PO perameters refer to 
the ionosphere's F layer • Using this models the refractive effects 
of the D and £ layers are not singled out» because they are quite 
small in comparison with those due to the F layer and are approxi- 
mately accounted for by allowing the electron density at the bottom 
edge of the F layer to be zero. 



(173) 
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Figure 24. Normalized 3-Parameter Model of Atmosphere 
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Computation of Ray Bending 

The major effect of ray bending due to refraction is that the 
observed elevation angle, E, of a vehicle in view of a station is 
greater that the actual or computed value of elevation, E^. This 
effect is clearly seen in Figure 25« The difference between the com- 
puted and observed elevation is the angle 6. 

In Figure 25 three vectors of importance are shown: 

p, the vector from the station to the vehicle, 

p'j a vector from the station toward the observed 

vehicle position (i»e», tangent to the ray 

path at the station), and 

■*■ -*■ 

Ap, a vector constructed perpendicular to p* 

The value of azimuth. A, and computed elevation, Ect are ob- 
tained from Equations (156) and (157) using the value of p, since 
these values do not incorporate refraction effects. From the geome- 
try of Figure 25, one can find the components of Ap in the east, 
north, and vertical station-oriented coordinate system. Thus: 



Ap - [q\ tan 6[T] 



sin Eq sin A 
sin E^ cos A 
- cos Eq 



where |p| is the magnitude of the vector p, and [T] is the 3^3 
matrix of Equation (15A) used to transform the coordinates of Ad 
from the station oriented coordinate system into the inertial system. 
From the figure it is also obvious that 



or 



p' + Ap 



p - \t\ tan 6[T] 



sin Ec sin A 
sin E^ cos A 
- cos E« 



In the program the vector p' is actually normalized to a unit vec- 
tor. The length of p ' is not important, however, since it is used 



(17M 



(175) 



(176) 
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Figure 25. Geometry of Ray Path Used to Obtain ^' 
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Figure 26. Geometry of Ray Path Used to Obtain dy 
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only to calculate observed angle-related measurements such as azimuth^ 
elevation, it-cosine^ m-cosine, X-angle^ Y-angle, and indirectly, hour 
angle and declination as described above* 

In Equation (176) all quantities are easily computed with the 
exception of 6^ the "error" in the elevation angle. The method of 
computing 6 used by the program is that due to Weisbrod and Anderson ■* . 
Only an outline of their analysis is presented here. 

First, consider the ray of Figure 26 entering an infinitesimal 
layer of thickness dr at an angle 6. At a boundary between two 
infinitesimal layers Snell's law holds. Thus; 

n cos 6 ■ constant (177) 

where n is the index of refraction. Differentiating Equation (177) 

with respect to path lengthy s^ yields Equation (178) below. Since 

n varies only with altitude or distance along an earth radius ^ r, 

one obtains: 

4^ cos 6 - n 4^ sin 6 - (178) 

as as 

dn dr ^ dS j o « f^ mn 

r— r— cos 6 - n -T" sin 6-0 (l/^^O 

dr ds ds 

but ~ - sin 6, (180) 

and from the definition of curvature one gets: 



(181) 



ds K 

where < is the radius of curvature of the ray. Thus, after substi- 
tution. Equation (179) becomes: 

i.i^^cos 6 . (182) 

K n dr 

From Figure 26» an infinitesimal length ds of the ray path is 
given by: 

ds - K dY - -rh (^^3) 

sin 6 



103 



thus I 

Combining Equations (184) and (182) r«ftult« in an •3cpr«ftaion for dyt 

dY - ^ ^ cct dr (185) 

AI8O9 from Figure 26, it is s««n that the dy** of all eleaantary 
layers are directly additive* Therefore, the angle Yj|^ due to the 
total bending between layers bounded by the heights Ta and rj^ is 
simply. the integral of Equation (185): 



Yjk - 



1 d[a 



^ ,^ cot 6 dr (IbO 

n or 



If a ray departs from a given layer (where r * r^ and n ■ n^c) ^t an 
angle of 6^, from Snail's law for spherical stratification, one has: 

nk ^k ^^* ^k ' constant (187) 

If the angle Yjic ^^* ^^ refractive bending is computed for only 
a very small layer of atmosphere between the heights ^ * ^j ^^^ 
r « r^, the following assumptions are justified: 

(1) — » constant I 

(2) the index of refraction n is close to unity, and 

(3) Ah ■ r;^ - rj « rj (Ah is an infinitesimal layer 
of height). 

Using these assumptions along with Equations (186) and (187), Weisbrod 

and Anderson *" -^ derive an expression for Yj;^ in terms of the angle 

6 and the refractivity N (see equation (168) for its definitian)^ 

for the details the reader is referred to their paper. Their expres«- 

sion for yjk through a small layer of atmosphere is: 

. yi^ ' 500 (ta!^j !\>n B,) (-^11^"'*^—) ("«> 

The total bending through the atmosphere Is the sua of the 
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individual contributions. Thus, 



= E 



h-i - \ 



.. 500 (tan P, - + tan p.) 
x=l '^i-l '^x 



(milliradians) 



(189) 



From Equation (187), the value of cos p is given as: 

i 



cos p. = 



n. . r. . cos p. . 
X-1 X"l x-1 

n. r. 

1 X 



n. . cos P . . 

1-1 X"l 

n. 

X 



Ah 



R + h. 
e X 



(190) 



and it follows: 



tan B . 
^x 



J. 



2 u 
cos p , 



COS P . 



(191) 



Equations (189), (190), and (191) determine the value of y 
where P . is the value of the angle of incidence of the ray at a 
height h. N. = N(h.) is the refractivity at height h. obtained 
from the atmospheric model, and Ah is the increment of height used 
in Equations (189) and (190). Note that h. = h._ + £^ and that 
h is the computed height of the vehicle obtainable from the program. 
At the station h = h^^, which is the height of the station, and P-. 
is set equal to the computed elevation, E . 

Once Y has been found, consideraticHi of the geometry in Figure 
27 leads to the evaluation of 6, the "error" in the elevation angle. 
Thus, we note from Figure 27: 



e = Y - (a - P) 



(192) 



Using Snell's law for spherical stratification and the application of 
the law of sines to triangle SOQ results in: 



COS P = "^ — COS p 



nr 







(193) 



r^ = pQ = r cos a 
Combining Equations (193) and (19A) yields: 



cos P =■ 



jj COS a = cos (a - (a 
105 



P)) 



(194) 



(195) 



y - (a+c) - I - (0+Y) 




Figure 27. Geometry of Ray Path Used to Obtain 6 
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Since a - P is a small angle, a small angle approximation and an 

approximation for "~r is used^ -^ in Equation (195) to obtain 

an expression for a - ^ in terms of refractivity at the station 

(Nq) and at the vehicle (N): 

a - 6 « (No - 10 10"^ cot P (196) 

Using Equations (192) and (196), the value of c is determined. 
Applying the law of sines to triangle SOP results in: 

ro COS (Bo - 5) - r cot ((jt + t) - M (197) 

Combining Equations (194) and (197) givet^ an expression for i; 

X itlu e tan a + (] ■ cos c) /,n-i\ 

tan 6 - 7^'" , ■ ^ ' ' ' ^ \ (197) 

•In c + cos c tan a - tan Pq 

Finally, using small angle approximations for 6 and f : 



^ £2 

£ tan a + ~ 
e -f tan a - tan Pq 



This value of 5 can now be used in Equation (176) to find the 
vector p', which determines all angle related observations includin; 
the effects of refraction. Due to the approximations used in the 
analysis, below elevations of about 5** the errors in the propagation 
corrections amount to a few percent. Because of this fdct* propapa- 
tion corrections due to refraction should be computed only when the 
elevation of the vehicle exceeds 5*. 

Refraction Effects on Range and Rang^e Rale 

Because of signal retardation caused by the refractive grr.ditrt 
of the atmosphere, the round-trip range computation of Equation (158) 
includes a correction term, Ap. The signal retardation caused by nn 
infinitesimal layer of thickness dr is given by: 



(199) 
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dT 






CSC B dr 



- (* - 1) 



CSC S 



dr 



-6 csc^^^ 



- N X 10 

c 

The effect on range between layers bounded by heights r - rj 
r - r^ is given by: 



(2(;'^) 



and 



Apjk 



^v ' 


•k 

* 


cdT " 

9 \ 


« 



N X 10"^ CSC 6 dr 



(2(.l) 



"^J 



^J 



Using methods similar to the computation of retractive bending, 
Weisbrod and Anderson ^ arrive at the expression for the rcfr*ctivt 
round-trip range correction, 2^p, due to the passage of the ray 
through t\ie entire atmosphere: 



2Ac 



2 X 10 



-6 



sin $1-1 + sin 6^ 



i-1 



+ 10"^ 



'^7f 



i-n |n^.) - N^Khj - \\_^ ^ 

^ sin 6i-i + sin 6^ 

i-m+i "• ^ 



C 



1. 



where 



^'i 

h, 



is the refractivity in the i 



th 



layer. 



f2 



is the height of the i^^ layer, 
is obtained from equation (190), 
is the transmitted (or up) frequency, 
is the received (or down) frequency. 

The indices in the first summation refer to layers in the troposphere, 
and those in the second refer to layers in the ionosphere, Th« values 
of U^ - N(hi) are obtained from the models of rcfractlvt index. 

Due to refractive bending, the range-rate measurement Of Equa- 
tion (159) also includes a correction term, Ap, This correction 
term arises because the program initially computes the component of 
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vehicle velocity along the vector p, i.e., along the straight line 
SP of Figure 27. The measured range-rate is the component of vehi- 
cle velocity along the tangent to the ray path at the vehicle, i.e., 
along the straight line TP of Figure 27. 

The velocity of the vehicle and the station is available from the 
program. If one denotes the component of vehicle velocity along the 
line SP by v^, and the component of vehicle velocity perpendicular 
to SP by Vy, then the range rate computation without the correc- 
tion for refraction (i.e., Ap - 0) is: 



Pc - ^1^ - vx (203) 

The measured component of velocity along the tangent to the ray path 
at the vehicle as seen from Figure 27 is 

pjjj « Vx cos (y - <S) + Vy sin (y - 6) (204) 

The correction to range-rate, Ap, is given by: 

Ap - Pn; - Pc (205) 

A factor should also be included to account for the round-trip fre- 
quency dependent effects in the ionosphere. After making a small angle 
approximation the final equation used for Ap , the one-way correc- 
tion to range-rate, is: 



(206) 



• 1 

Ap- ■2 


tl 


Vy (Y - -5) 



where 



f2 



is the component of vehicle-station velocity perpendicular 
to the vector p, and lying in the plane which includes 
the center of the Earth, 

is the transmitted frequency, 

is the received frequency. 
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SECTION III. 

USER'S GUIDE 
INTRODUCTION 

The information in this section is intended as a guide for those 
desiring to use SPACE-A. SPACE-A generates a vehicle ephemeris and 
associated observations from ground stations and has four modes of 
operation. They are: 

(1) Computation of vehicle ephemeris only; 

(2) Generation of obseirvations from ground stations; 

(3) Same as (2) with the writing of the observations onto 
tape in a format suitable for processing by SPACE-B; 

(4) Same as (2) with the following included as standard 
outputs: 

(a) the acquisition time, 

(b) total time in sight, and 

(c) the option of the printing time of polar and 
meridian crossings. 

DESCRIPTION OF INPUT DATA 

All input data is read at the beginning of each run. The format 
and arrangement of the data is designed to make the data prepara- 
tion as convenient as possible. It is seldom necessary to rend nil 
the data; in fact, most runs require a minimum of cards. Many vari- 
ables have standard values which are preset by the program and should 
be satisfactory for most cases. 

The input data is divided into three parts. The first consists 
of a single card, called the basic input card, which must be included 
in each run. The second and third parts are designated group I and 
group II inputs, respectively. Each group is divided into several 
sections, each of which must be preceded by a card containing the 
section number, in integer format, in columns four and five. Any 
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section may be omitted from the input deck. Parts of sections, how- 
ever, must not be omitted. See Figure 28 for the setup of a sample 
input deck. 

The basic input card, which precedes each case, contains four 
variables. The first pertain$ to the standard values which may be 
set in place of the group II data and will be explained later. The 
second, MDE, indicates the function to be performed (i.e., trajectory 
computation, observables, etc.). The third indicates the desired pre- 
cision and determines the set of standard values used for integration. 
The last specified Encke or Cowell integration. 

The group I inputs consists of several sections. Any of these 
sections may be omitted from the input deck if the variables contained 
in the section are not needed for the run or in the case of stacked 
cases, the values from the last case are to be used again. The des- 
cription of each variable given in the data summary (see Table IV) is 
sufficient to understand most of the inputs. Therefore, only a brief 
discussion of the variables given in Table IV is presented here. 

In the group I input data, if the variable KLM2 of Section 2 is 
or 1, the initial conditions are assumed to be referred to the base 
date system or true system of date, respectively. 

Sections 7 and 9 deal with observations from observing systems 
on vehicles and with powered flight and have been omitted since they 
are not presently used. The variables in Section 8 control the out- 
put and are discussed in DESCRIPTION OF OUTPUT INFORMATION, 

The end of the group I input data is indicated by a card con- 
taining a 10 in columns four and five. This card must always be in- 
cluded in the data deck. After reading this card the program advances 
to the group II inputs. 

Before the program reads the group II input data, the val»e of 
the first variable contained on the basic input card, KSTDRD, is 
tested. If this variable is negative, standard valuts of the vari- 
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ablcs In the group II Input data ara computad from atorad Information 
and tha program than baglna to raad the group II Inputa to augment or 
replace any of the atandard values* Table V glvea a Hat of the atan- 
dard values uaed by tha program. If KSTDRD la not negative , none of 
the standard values is set and all valuaa should be raad In unleas 
values from a preceding stacked caae are used* Tha and of group II 
inputs Is Indicated by reading a card with a 20 In Integer format In 
columns four and five. 
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Table IV. 

Summary of Input Data 

Basic Input Card 



COLUMN 



2-5 



6-10 



NAME 
KSTDRD 

MDE 



FORMAT 



15 



15 



11-15 



PRECIS 



F5.0 



16-20 



CEPID 



F5.0 



DESCRIPTION 

< Will want at least some 
standard inputs from 
group II 

^ No standard. All values 
will be read in. 

- 1 Trajectory computation. 

■ 2 Observable computations 

without simulated data. 

■ 3 Observable computations 

with simulated data. 

Prediction mode 

Low precision level. 

Intermediate precision 
level. 

High precision level 

1 Encke integration 

1 Cowell integration 



- A 

- 1 

- 2 

- 3 

- 0, 



115 



Group I Inputs 



SECTION CARD 


COLUMNS 


NAME 


FORMAT 


DESCRIPTION 


1 1 


I'll 


ITITLE(1-12) 


12A6 


Title 


2 


2-5 


NYEARP 


15 


Year of initial 
conditions 




6-10 


DAYS 


F5.0 


Day of initial 
conditions. 




11-15 


HRS 


F5.0 


Hour of initial 
day. 




16-20 


HMIN 


F5.0 


Minute of initial 
hour. 




21T40 


SEC 


E20.16 


Seconds of initial 
minute. 


3 


1-2A 


TMAX 


E2A.16 


Time of run in 
hours. 


2 1 


2-5 


KLM 


15 


Indicator for units 
of position and 
velocity vector. 

- 1 KM, KM/SEC 

- 2 ER,ER/HR 

- 3 FT, FT/SEC 

- A MI,MI/HR 

- 5 NM,NM/HR 

- 6 NM, FT/SEC 

- 7 AU,AU/HR 




6-10 


KLMl 


15 


Indicator for coor- 
dinate system of 
input vectors. 



" 1 Cartesian coord. 
WE - 

- 2 Cartesian coord. 

Compute WE 

- 3 Geodetic long, 

lat, alt, 
VeiVn.Vh; WE • 

- A Geodetic long, 

lat, alt, 
VeiVnpVh; 
Compute WE 
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Group I Inputs 
(continued) 
SECTION CARD COL UMNS NAME FORMAT 



11-15 KLM2 



15 



16-20 KSNAP 



15 



21-25 KLM3 



15 



DESCRIPTION 



- 5 


Geodetic long. 




laty alt. 




V|, Y, az; 




WE - 


- 6 


Geodetic long. 




laty alt. 




|V|. Y. az; 




compute WE 


- 7 


Geocentric ra. 




declf alty 




VratVdtVh; 




WE - 


- 8 


Geocentric ra. 




decly alty 




Vra,Vd,Vh; 




compute WE 


- 9 


Geocentric ra. 




decly alty 




|V , 6, az; 




WE - 


- 10 


geocentric ra. 




decly alty 




V , 6, az; 




compute WE 


Indicator for nuta- 


tion 


and precession 


of input vector. 


- 


no 


- 1 


yes 


Indicator for nuta- 


tion 


and precession 


of vectors during 


run. 


. 


- 


no 


- 1 


yes 


Indicator for libra- 


tion 


of input vectors. 


- 


no 


- 1 


yes 
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Group I Inputs 



SECTION CARD 



3 
4 





(continued) 






COLUMNS 


NAME 
KLIBR 


FORMAT 
15 


DESCRIPTION 


26-30 


Indicator for libra' 








tlon of vectors 








during run 








- no 








- 1 yes 


31-35 


MRREF 


15 


Indicator for Ini- 
tial reference body 

- 1 Earth 

- 2 Sun 

- 3 Moon 

- 4 Mars 



1-24 

25-48 

49-72 
1-24 

25-48 



6-10 



RCIN(l) E24.16 

(2) E24.16 

(3) E24.16 
RCIN(l) E24.16 



(2) 



E24.16 



KSPLT 



15 



- 5 Venus 

- 6 Jupiter 

- 7 Saturn 

Initial position 
vector 

See KLM and KLMl 
for units and type. 



Initial velocity 
vector. 

See KLM and KLMl 
for units and type. 



49-72 


(3) 


E24.16 




2-5 


PASFX 


F5,0 


Total number of 
passes. 


2-5 


KS2BY 


15 


Indicator for two- 
body Integration 

- no 

- 1 yes 



Indicator for in- 
clusion of planetary 
perturbations 

- no 
• 1 yes 
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SECTION CARD 





Group 1 Inputs 






(continued) 




COLUMNS 


NAME FORMAT 
KSOBL 15 


DESCRIPTION 


11-15 


Indicator for Inclu 






sion of oblateness 






perturbation* 






- no 






- 1 yes 



16-20 



KSDRG 



21-25 



KSRAP 



26-30 



KSDRGM 



31-35 



KSDRGV 



36-40 



KSMNOB 



41-A5 



KRF 



15 Indicator for inclu- 
sion of Earth drag 
perturbation, 

- no 

■ 1 yes 

15 Indicator for inclu- 
sion of radiation 
pressure perturba- 
tion* 

- no 

■ 1 yes 

15 Indicator for inclu- 
sion of Mars drag 
perturbation, 

- no 

- 1 yes 

15 Indicator for inclu- 
sion of Venus drag 
perturbation, 

- no 

■ 1 yes 

15 Indicator for inclu- 
sion of Moon oblate- 
ness perturbation. 

- no 

■ 1 yes 

15 Indicator for inclu- 
sion of refraction 
effects. 

- no 

■ 1 yes 
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Group I Inputs 



SECTION CARD 





(continued) 






COLUMNS 


NAME 


FORMAT 


DESCRIPTION 


46-50 


KECLPS 


15 


Indicator for inclu- 
sion of eclipse in- 
formation print. 

- no 

- 1 yes 


2-5 


MAXSTA 


15 


Total number of 
stations used in 
run. 


2-5 


K 


15 


Station number , 


10-15 


STANM(L) 


A6 


Station name. 


20-30 


TYPE(L) 


4X111 


A(1)+1.E2*A(2) 



+l.EA*A(3)+i.E6*A(A) 

+1,E8*K where 

A - observation 
types used by 
station K in 
ascending order, 

L - packed station 
number 



1-24 


STALT(K) 


E24.16 


Latitude of station 
K. 


25-48 


STALN (K) 


E24.16 


Longitude of station 
K. 


49-72 


STAHT (K) 


E24.16 


Altitude of station 
K. 


1-36 


RRATE(1,L) 


E36.16^ 


[^Repetition rates 
(hrs). 


37-72 


RRATE(2,L) 


E36.16 


\ \ For each observation 


1-36 


RRATE(3,L) 


E36.16 


LType 


37-72 


RRATE(4,L) 


E36.16J 




1-18 
19-36 
37-54 
55-72 


TDELAY(1,L) 
TDELAY(2,L) 
TDELAY(3,L) 
TDELAY(4,L) 


E18.8 
E18.8 
E18.8 
E18.8 


r Times in hrs before 
j which each observa- 
S tion is not to be 

V computed. 
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SECTION CARD COLUMNS 

7 1-18 

19-36 

8 1-2A 
25-48 
49-72 

9 1-24 
25-48 

49-72 



10 



11 



12 



1-24 

25-48 

49-72 

1-24 

25-48 

49-72 

1-24 

25-48 

49-72 



Group I Inputs 
(continued) 

NAME FORMAT 

FUP(K) E18.8 

FDOWN(K) E18.8 

STA0R(NCDST+1) E24.16 

STA0R(NCDST+2) E24 . 16 

STA0R(NCDST+3) E24.16 

STA0R(NCDST+4) E24.16 

STA0R(NCDST+5) E24.16 

STA0R(NCDST+6) E24 . 16 

STA0R(NCDST+7) E24.16 

STA0R(NCDST+8) E24 . 16 

STA0R(NCDST+9) E24 . 16 

STAOR (NCDST+10) E24 . 16 

STAOR (NCDST+11) E24 . 16 

STA0R(NCDST+12) E24.16 

STAOR (NCDST+13) E24 . 16 

STAOR (NCDST+14) E24a6 

STA0R(NCDST+15) E24 . 16 



DESCRIPTION 

Station transmit 
freq. (MHz). 

Station receive 
freq, (Mllz). 

AEE station 
rotation angle. 

AEV station 
rotation angle. 

AEN station 
rotation angle. 



AU 
AV 
AW 



r. 



station lo- 
cation 

errors 
caused by 

geodetic 
net error. 



NO; refractiv- 
vity at sea 
level, 

H; troposphere 
scale fact, 

PC; max. elec- 
tron density. 

HO; loser limit 
of ionosphere. 

HM; Ht of PO 
(KM) 

- open - 

AT for timing. 

Bias added for 
obser.A(l) 

Bias added for 
obser.A(2) . 
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Group I Inputs 












(continued) 






SECTION 


CARD 
13 


COLUMNS 
1-24 


NAME 
STA0R(NCDST+16) 


FORMAT 
E24.16 


DESCRIPTION 




Bias added for 
obser,A(3) 






25-48 


STA0R(NCDST+17) 


E24.16 


Bias added for 
obser.A(4) 






REPEAT 


CARDS 2-13 FOR EACH 


STATION 




6 


1 


1-18 


DAREAS 


E18.8 


Effective sur- 
face area of 
vehicle per- 
taining to 
drag (ft2). 



9 

10 



19-36 



PAREAS 



37-54 

OPEN 
1 



SPADD(6) 



E18.8 Effective sur- 
face area per- 
taining to 
radiation pres- 
sure (ft^). 

E18.8 Mass of vehicle 
(lb-masses) 



2-5 


lUNIT 


15 


Indicator for 
output units 
(see KLM) 


6-55 


IPSEC(I),I-1.10 


1015 


Indicator for 
suppression of 
each of 10 print 
sections. 


1-18 


DTP I 


E18.8 


Print portion 


19-36 


DTSUP 


E18.8 


(hrs) and sup- 

^4__ 



37-54 



OPEN 



PRATE 



E18.8 



(hrs), ot total 
print period. 

Interval with- 
in DTPI at which 
to print. 



END OF GROUP I INPUTS 
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SECTION 


CARD 


COLUMNS 


NAME 


FORMAT 


DESCRIPTION 


X 


1-7 


1-72 


(DT3(I»J), 












1-1,3) J-1,7 


3E2A.12 


Integration in- 
tervals for each 
of 7 working 
bodies for near, 
medium and far 
reference. 


2 


1 


1-72 


Rl(l-3) 


3E2A,12 


Rl and R2 are 




2 


1-72 


Rl(A-6) 


3E2A.12 


distances in 
E.R. for each of 




3 


1-72 


Rl(7),R2(l-2) 


3E2A.12 


7 working bodies 




A 


1-72 


R2(3-5) 


3E2A.12 


for switching 
from near to 




5 


1-A8 


R2(6-7) 


3E2A.12 


medium and me- 
dium to far in- 
tegration inter- 
vals. 


3 


1 


1-2A 


RTl 


E2A.12 


Values used as 






25-A8 


RT2 


E2A.12 


tolerances in 



1-18 



19-36 



DHl 



DH2 



criteria. 

E18 • 8 Troposphere 
Integration 
step size (KM) » 

E18.8 Ionosphere in- 
tegration step 
size (KM). 



37-5A 


H2 


E18, 


,8 


Upper limit 
troposphere 
(KM). 


55-72 


HA 


E18, 


,8 


Upper limit 
ionosphere (KM) . 


2-5 


KOBLAT 


15 




Number of ob- 
lateness coeffi- 
cient terms. 


2-5 


N 


15 




N index. 


6-10 


M 


15 




M index. 


1-2A 


SORCl 


E2A.12 


Value of C 
coefficient. 
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Group II Inputs 







(continued) 






SECTION CARD 


COLUMNS 


NAME 


FORMAT 


DESCRIPTION 




25-48 


S0RC2 


E2Aa2 


Value of s 
coefficient. 


REPEAT 


CARDS 2-3 


UNTIL KOBL^T VALUES HAVE BEEN 


READ IN, 


6 1 


2-5 


MBMAX 


15 


Number of working 
bodies. 




6-10 


KWBMU(I) 
I*l, MBMAX 


15 


Indices of 
working bodies. 


2 


1-72 


TPMAT8(1) 
I-l, MBMAX 


3E2A.12 


Gravitational 
constants of 
working bodies. 


REPEAT CARD 


2 FOR EACH VALUE 


OVER 3N NEEDED. 


7 1 


1-2A 


DYN(48) 


E2A.12 


Solar flux in 
10-22 w/m2-Hz 
at 10.7 cm. 




25-48 


DYN(49) 


E2A.12 


Open 


8 1 


1-72 


DYN(31-53) 


3E2Aa2 


Coefficients 
for lunar ob- 
lateness 
(kg • km2). 


9 1 


1-24 


COMB(l) 


E2Aa2 


Velocity of 
light. 


10 1 


1-72 


PRKT3(l-3) 


3E2Aa2 


Print intervals 
(hrs) for near 
medium and far 
reference. 



11 



12 



1-2A 



1-18 



EMlN E2A,12 Minimum value 

of elevation 
angle (radians) . 

RTO E18,8 Ratio of 

Nordsieck inte- 
gration interval 
to that in Runge- 
Kutta. 
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Group II Inputs 









(continued) 






SECTION 


CARD 


COLUMNS 


NAME 


FORMAT 


DESCKIPTION 


13 


1 


1-24 


DSPL 


E24.12 


Special inte- 
gration inter- 
val in A4 mode 
to obtain ac- 
quisition time 
(hrs). 


14 


1 


1-72 


RATEC( 1-3,1) 


3E24.12 


Rotation vector 
used in Mars 
drag computa- 
tions. 




2 


1-72 


RATEC(l-3,2) 


3E24.12 


Rotation vector 
used in Venus 
drag computa- 
tions. 


15 


1-10 


1-72 


XMACH(I) 
1-1,40 


4E18.8 


Mach number 
table. 




11-20 


1-72 


CDT(I) 
1-1,40 


4E18.8 


Drag coefficient 
table 


16 


OPEN 










17 


OPEN 










18 


OPEN 










19 


OPEN 










20 


END OF 


GROUP II 


INPUTS 
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Table V. 
Standard Values for Group II Inputs 

Section 1, 

DT3 ARRAY (3^7 array of integration intervals) 



1-1,3, A, 5, 6, 7 ,1953125 x lO^^ 

1-1,3, A, 5, 6, 7 ,15625 x 10"^ 

1-1,3, A, 5, 6, 7 .125 



DT3(1,I) 

DT3(2,I) 

DT3(3,I) 

DT3(1,2) - A. 

DT3(2,2) - 6. 

DT3(3,2) - 10, 

Note: The above values are always used for Encke integration and for 
Cowell when low precision is specified* When Cowell and inter- 
mediate precision is specified, each value in the DT3 array is 
divided by 2; when Cowell and high precision is specified each 
value is divided by A. The precision desired is controlled by 
the variable CEPID on the Basic Input Card. 

Section 2. 

Rl ARRAY (1x7 array of reference body change 

criteria) 

R1(I) - A • RADII(I) 

R2(I) - 10 • RADII (I) 
where RADII (I) is the radius of the Ith working body (Earth, Sun, 
Moon, Venus, Mars, Jupiter, Saturn). 

Section 3. 

RTl - 1. X 10-^ 

RT2 - 1. X 10"^ 
If low precision is specified 

RTl - 1 X 10-** 

RT2 - 1 X 10-** 
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Table V. 
(continued) 

Section 4. 

DHl - 1. 
DH2 - 10. 

H2 - 40. 

H4 - 2,000. 

Section 5. 

KOBLAT - 3 

C20 - - 1082.3 X 10~^ S20 - 
C30 - 2.3 X 10"^ S30 - 
Ci»o - 1.8 X 10"^ Si»o - 

Section 6. 

KBMAX - 7 

KWBMU(I) - I; I - 1,2,3,4,5,6,7 

TPMAT8(1) - 19.909416 

(2) - 6629965.8 

(3) - . .24478289 

(4) - 16.1983009 

(5) - 2.14364682 

(6) - 6338.16258 

(7) - 1897.36734 

Section 7. 

DYN(48) - 200. 

Section 8. 

DYN(51) - .88746 " lO^^ kg • ko^ 

DYN(52) - .88764 « lO^' kg • kn^ 

DYN(53) - .88807 * lO^^ ^g • kr^ 
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Table V. 
(continued) 

Section 9. 

COMB(l) - 169210. 5800A 

Section 10. 

PRNT3(1) - .25 
PRNT3(2) - .5 
PRNT3(3) - 1. 

Section 11. 

EMIN - (5. degrees; expressed in radians) 

Section 12. 
RTO - 3 

Section 13. 

DSPL - 1. X 10'3 

Section 14. 

RATEV(1,1) - 0.0 

RATEV(2,1) - 0.0 

RATEV(3,1) - 0.255175469 

RATEV(1,2) - 0.0 

RATEV(2,2) - 0.0 

RATEV(3,2) - 0.0 
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Table V, 
(concluded) 



Section 


15. 






MACH TABLE 


ARRAY 


(XMACH(I), I 


- .5 


1.5 


25. 


96. 


0.0 


1.85 


30. 


97. 


.5 


2.0 


50. 


98. 


.75 


■2.7 


70. 


99. 


.85 


3.4 


90. 


100. 


.9 


4.2 


91. 


101. 


1.04 


5.6 


92. 


102. 


1.1 


6.75 


93. 


103. 


1.2 


8.5 


94. 


104. 


1.3 


15. 


95. 


105. 


DRAG COEFFICIENT TABLE (CDT(I 


.8 


1.4 


1.14 


1.14 


.8 


1.38 


1.14 


1.14 


.82 


1.36 


1.14 


1.14 


.92 


1.285 


1.14 


1.14 


.99 


1.23 


1.14 


1.14 


l.OA 


1.19 


1.14 


1.14 


1.18 


1.16 


1.14 


1.14 


1.21 


1.155- 


1.14 


1.14 


1.26 


1.15 


1.14 


1.14 


1.3 


1.45 


1.14 


1.14 



1,40) 
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DESCRIPTION OF OUTPUT INFORMATION 

The input quantities which determine the form of the output of 
the program are given in Section 8 of the Group I inputs in Table IV. 
The choice of output units are determined by lUNIT. The choices of 
units with corresponding values of lUNIT are given below. 



Value of lUNIT 


OUTPUT UNITS 


1 


KM, KM/SEC 


2 


ER,ER/HR 


3 


FT, FT/SEC 


4 


MI,MI/HR 


5 


NM,NM/HR 


6 


MM, FT/SEC 


7 


AU,AU/HR 



Time may be divided into periods of print and suppression of 
print by setting the input variables DTPI and DTSUP. A period of 
printing (DTPI) is executed followed by a period of suppression (DTSUP) 
and the cycle is then repeated. 

The variable PRATE, is used to determine the print interval with- 
in DTPI. Since conditions occur where the print routine is entered 
more than once with the same time, printing is normally only executed 
upon the first entry. However, by using a negative PRATE, printing 
will occur each time the print subroutine is entered during the print 
interval determined by DTPI. 

The user has the option of choosing the information he would like 
printed by setting IPSEC(I) « 1, where I is the output section number 
given in Table VI, that contains certain information. By reading 
a zero for IPSEC(I), the printing of the I^h section is surpressed. 
There are ten sections of output information as shown in Table VI. See 
Figure 29 for a sample listing with various sections labeled. 
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Table VI. 
Ten Sections of Output Information 

Section 1. 

Program time in hours since launch; 

Current program date in hours, minutes, seconds, 
day and year; 

Planetary body which is the current reference. 

Section 2. 

Station number; 

Observation types which station handles. 

Section 3. 

Components and magnitudes of position and velocity 
vectors in units specified by the user. 

Section 4. 

Components and magnitudes of perturbations in 
position, velocity and acceleration in units 
specified by the user. 

Section 5. 

Components and magnitudes of vectors between 
vehicle and each planetary body in specified 
units. 

Section 6. 

Components and magnitudes of vectors between 
the Earth and each of the other planetary bodies 
in specified units. 

Section 7. 

Components and magnitudes of vectors between 
the Sun and each of the planetary bodies 
in specified units. 
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Table VI. 
(concluded) 

Section 8. 

Right ascension; 

Declination; 

Greenwich hour angle; 

Longitude and latitude of Earth subsatellite point; 

Geocentric altitude, azimuth and elevation; 

Geodetic azimuth and elevation. 

Section 9, 

True, eccentric and mean anomalies; 

Semi-major axis; 

Eccentricity; 

Inclination; 

Argument of perigee; 

Perigee passage time; 

Right ascension of ascending node; 

Period; 

Mean motion; 

Perigee and apogee heights; 

Unit vector to perigee; 

Unit angular momentum vector. 

Section 10, 

Station name; 
Program time in hours; 
Station location; 
Observation values. 
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PROGRAM DECK SETUP 

This section describes the deck setup for running the A-tnode of 
SPACE on the 7030 computer. The program is executed under the control 
of MCP (Master Control Program). The binary cards containing the com- 
piled routines are kept on a magnetic tape at the MITRE 7030 computer 
facility and only a few control cards and the data deck are needed by 
the user. 

Those cards shown in Figure 30 that have a B or T punched in 
column 1 are control cards for the MCP, They are described in greater 
detail in the MITRE 7030 Facility Manual. 

JOB Card 

The job card is always the first card of an input deck. This 
card contains necessary information for accounting and no job can be 
run without it. The fields on the JOB card are variable in length 
and separated by commas. 



Format 

1 
B 



10 

JOB, 'JQBID' /NAME'/ PROJECT '/DEPT*/MAXTIME'/DEST' 



Field Name Description 

B Punched in column 1 

JOB, Punched in columns 10-13 

'JOBID' 1 to 7 character job identi- 

fication 

'NAME' ' Nam« of programmer (up to 13 

characters) 

'PROJECT' Project number (3 to 9 

characters) 

'DEPT' Department number (3 

characters) 

'MAXTIME' Maximum running time in 

minutes (up to 3 characters) 

'DEST' Destination of output 

(up to 6 characters) 
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Example: 

1 

B 



TYPE Card 



10 

JOB, SPACE, BOND J, 007, D84, 13, AlAl 



This card defineB the type of job (GO), and the proper monitor 
system (FORTRAN) to the MCP to execute this program. 



Format: 

1 

B 



10 

TYPE, GO, FORTRAN 



ADDIO Card 

This card contains the library number of the program tape that 
contains the binary decks of the SPACE program* Since modifications 
to the program may be made in the future, the user should periodically 
check the current tape number with Department D-8A. 



Format; 

1 
BLIBl 



10 

ADDIO, PLBXXXX 



Field Name 

BLIBl 

ADDIO, PLBXXXX 



Description 

Punched in columns 1-5 

Beginning in column 10, 
the XXXX indicate the 
field where the tape 
number is to be punched* 



SUBTYPE, BIN Card 



This card tells the MCP that the cards that follow are binary* 



Format : 

1 
T 



10 
SUBTYPE, BIN 



TMAIN Card 

This card results in the construction of a dummy main program 
which calls the executive routine of SPACE (EXECA). EXECA is compiled 
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as a subroutine and the binary cards are contained on the library tap«« 
This card Is the equivalent to the Fortran routine. 

EXECA 







CALL 




RETURN 




END 


Format: 




1 


10 




T 


MAIN, EXECA 




Processed FIOD Deck 





Since the programmer cannot choose the Input-output channels to 
be used (they are assigned by MCP at execution time), the FIOD deck 
specifies the I/O channel requirements of the program. The processed 
FIOD deck Is the result of compiling the following subprogram: 

1 10 

T SUBTYPE, FIOD 

B2 IOD,$READER 

B3 IOD,$PRINTER 

B9 lOD, TAPE,,,,,, SAVE 

B8 lOD, TAPE,,,,,, SAVE 

B REEL,PLBDL680 
END 
Logical units 2 and 3 are the system Input and output units res- 
pectively. Logical unit 9 Is used for the output of the observation 
tape generated when that option Is exercised. Logical unit 8 Is used 
for the Input of the ephemerls tape containing the positions of the 
Moon and planets. The reel card specifies that the ephemerls tape, 
whose label Is PLBL680, Is to be mounted* 

SUBTYPE. DATA Card 

This card tells the MCP that the following cards contain the 
data for the program. 
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Format: 

1 
T 



10 
SUBTYPE, BIN 



Input Data Deck 

This Is the data deck described in the Description of Input Data. 



(INPUT DATA DECK) 



I 



SUBTYPE, DATA 



(PROCESSED FIOD DECK) 



MAIN, EXECA 



T 



SUBTYPE, BIN 

ADDIO, PLBXXXX 



B TYPE, GO, FORTRAN 





Figure 30. Deck Setup 
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SECTION IV. 
PROGRAMMING 

GENERAL FLOW OF PROGRAM 

The general flow of the program is shown in Figure 31« The user 
does not need to specify the number of cases to be run, as the program 
returns to read the data for the next case until all data cards are 
exhausted* 

The program structure and linkage of the program between basic 
routines is shown in Figure 32. It is to be noted, however, that the 
figure does not indicate every linkage but only the major ones. For 
instance, in some cases it is necessary for OBSERA (which usually only 
computes observations) to call the integration package. 

The routines listed under the heading of General Routines, in 
Figure 32 are used by many of the major routines. The linkage of 
these routines is not included as this would defeat the purpose of 
the figure, which is to indicate the position of the major routines 
in the overall structure of the program. 

DESCRIPTION OF SUBROUTINES 

The following part of the document gives a brief description of 
the subroutines used by SPACE-A. Each subroutine description gives 
the purpose and a brief description of the subroutine, as well as the 
other subroutines which call and are called by the given subroutine. 
For understanding the details of each subroutine one ought to refer 
to the pertinent portions of Section II and Section III dealing with 
the function of that subroutine as well as referring to the computer 
listing. 

It should be noted that, as of the writing of this document, 
certain subroutines have not been tested. These non-operational sub- 
routines are: LUNOBL, MVDRAG, OBD, PFINIT, PFLGHT, STACUL. These 
subroutines deal with lunar oblateness gravitational perturbations. 
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RETURN TO 
FIND NEXT 
ACTIVITY 
TIME 




INITIALIZE AND READ 
INPUT DATA 



NO 



1 



PICK NEXT 
ACTIVITY TIME 



PERFORM INTEGRATION 
(ENCKE or COWELL) 



COMPUTE OBSERVATIONS 
(IF REQUIRED) 



PRINT DESIRED 
INFORMATION 




YES 



RETURN 
TO READ 
DATA FOR 
NEXT CASE 



Figure 31 . General Flow of Program 
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EXECA 
(EXECUTIVE ROUTINE) 



T 



INPUT A 
(READS DATA) 



I 



MAIN A 

(CONTROLS MAIN 
FLOW OF PROGRAM) 



XFORM 

(READS a CONVERTS 
INITIAL CONDITIONS) 



I 



TIMNGA 

(FIND NEXT 
ACTIVITY TIM^ 



I 



PFINIT 

(INITIALIZE _ 
POWERED FLIGHTl 



STACUL 

(CHECK FOR 
OCCULATION) 



OBSER A 

(COMPUTES 
OBSERVATIONS) 



PRINT A 

(PRINT SELECTED 
INFORMATION) 



COW ELL 




INTEGRATION 
ROUTINES 




OBD (ON-BOARD OBSERVATIONS) 



• ATIM (MODE 4 COMPUTATION) 



• STAPOS (STATION VECTOR) 



ENCKEX * MODELA (INDEX OF REFRACTION) 



CITGRA ( CONTROLS INTEGRATION) 
CCHREF ( CHANGE REFERENCE BODY) 
CRSTRE ( SAVE /RESTORE VECTORS) 
CINT (INTEGRATION) 



GENERAL ROUTINES i 

DDOT (COMPUTE A-"b) | 

DMTML (COMPUTE H"[b]/H-[b] )' 



DOMUD (CHECK FOR ERROR) I 

I 
EPHEM (READ EPHEMERIS TAPE^i 

FIX (UNPACK WORD) 

NUTPRE (COMPUTE NUTAT/ 
PREC.) 

SERVCE (COMPUTE AX B, 

ICI) 



EITGRA (CONTROLS INTEGRATION) 
ECHREF ( CHANGE REFERENCE BODY) 
ERSTRE (SAVE/RESTORE VECTORS) 
EINT (INTEGRATION) 

KEPLER (N0MINAL/2-B0DY ORBIT) 
RECT (RECTIFICATION OF ORBIT) 



I 



DERIV 
(COMPUTES ACCELERATIONS) 



I 



OBLDRG 

(OBLATENESS 
AND DRAG) 



(Vl 
(D 



DENSTY 

(COMPUTES 
DENSITY 
OF AIR) 



I 



LUNOBL 

(LUNAR 
OBLATENESS) 



ECLIPSE 

(CHECK 8k PRINT 
ECL. INFO.) 



I 



MVDRAG 
(MARS/VENUS; 



PFLGHT 

(POWERED 
FLIGHT) 



Figure 32. Program Structure 
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Mars and Venus drag model, on-board observation*^ parturbatlona due 
to vehicle thrust, and observations dealing with occultation. 

EXECA 

Purpose ; 

This is the executive routine for SPACE-A. 

Description ; 

The program alternates between calling INPUTA and MAINA until 
all cases have been exhausted. 

Subroutines Called : INPUTA, MAINA 

ATIM 

Purpose: 

This subroutine is used in the visibility computation mode (mode 
4)« It determines the time a vehicle comes within sight of a given 
ground station or the time of polar baseline crossing (the time that 
the vehicle crosses the east~west vertical plane through the zenith) 
and meridian crossing* 

Description : 

When the vehicle is in sight, the i and m direction cosines 
are tested and used to interpolate for the times of polar and meri- 
dian crossings. When the acquisition time is to be determined, the 
integration routines CINT and EINT or the two-body routine KEPLER is 
employed to predict ahead until the vehicle comes into view. 

Called By : OBSERA 

Subroutines Called : CINT, CRSTRE, EINT, ERSTRE, KEPLER 

CCHKEF 

Purpose ; 

This subroutine tests criteria for changing the reference bo<iy 
when the Cowell integrator is used. Even though a change of reference 
body in the Cowell method is not necessary, it is utilised in the 
program in the same manner as in the Enc^ integrator. 
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Descrlptlon t 

Criteria for changing the reference body based on regions of 
Influence are tested* 
Called By : CITGRA 
Subroutines Called : DDOT, EPHEM^ SERVCE 

CINT (EINT) 

Purpose : 

This subroutine perforins the integration when the Cowell formu- 
lation is employed* 

Description : 

The Gill modification of Runge->Kutta is used for short-^term inte- 
gration and as a starting procedure for the Nordsleck long->term inte- 
gration* 

Called By : ATIM» CITGRA, MAINA, OBSERA 

Subroutines Called : CCHREF, CINT, CRSTRE 

CITGRA 

Purpose : 

This subroutine serves as the sub-main program governing calls 
to the integration routines in the Covell method* 

Description : 

If not in powered flight, the program checks for a change of re- 
ference bodies* The delta of integration is determined by the dis- 
tance of the vehicle from tht refetence body and integration is per- 
formed up to the next activity time by calling the integration routine* 

Called By : MAINA 

Subroutines Called : CCHREF, CiNt, CRSTRE 

CRSTRE (ICR) 
Purpose : 
This subroutine saves or restbrcs time, position, and velocity 



145 



and acceleration of the vehicle at any designated time for the Cowell 
integrator. 

Description ; 

Depending on the value of the argument (IRC) the information is 
saved (IRC <^ 1) or restored (IRC > 1). 

Called By ; ATIM, CITGRA, MAINA, OBSERA 

DDOT(A.B) 

Purpose: 

This function computes the dot product of two vectors. 

Description ; 

The program uses the formula 

DDOT - A(l)* B(l) + A(2)* B(2) + A(3)* B(3) 

Called By ; CCHREF, DERIV, ECHREF, LUNOBL, MVDRAG, OBD, 
OBLDRG, OBSERA, PRINTA, RECT, STACUL 

DENSTY(X1, X2, X3, DENP, DENEQ) 

Purpose : 

This subroutine evaluates the density of air in the upper 
atmosphere. 

Description ; 

Tables of the logarithm of p • C^^ obtained from the Harris- 
Priester data are stored in this routine. Interpolation is performed 
to calculate log(p • C^y) (DENEQ) at the equator as a function of the 
altitude of the vehicle (XI), solar flux (X2) , and solar time (X3) , 
Interpolation is also performed to calculate log(p • Cav)(DENP) at 
the poles, as a function of the altitude and solar flux. 

Called By: OBLDRG 



DERIV 



Purpose: 

This subroutine evaluates the acceleration terms for the Encke 
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and Covell Integrators • 

Description ; 

This subroutine computes the planetary perturbations and the solar 
radiation pressure perturbations as presented in the theory section of 
this manual. Earth oblateness perturbations and drag are computed by 
calling OBLDRG* Provisions are included for computing powered flight 
accelerations. 

Called By ; CINT> EINT 

Subroutines Called ; DDOT^ DOMUD^ ECLIPS^ EPHEM^ KEPLER » 

LUNOBL, MVDRAG^ OBLDRG, PFLGHT, SERVCE 

DMT^tL (A^ By C, I, J, Ky L» M^ N^ lAC, JAB^ KBC, IFLAG) 

Purpose; 

This subroutine multiplies two matrices of any size (up to 26 x 
26), or two matrices in which the second is transposed* 

Description ; 

The matrices A| B| and C are dimensioned I x J, K x L, and 
L X M respectively. 

When IFLAG " 0| 1, multiplication is done by rows of A so 
that A can be overwritten if desired (i.e., the argument C in the 
calling sequence may actually be the same as A) • 

When IFLAG " 2| 3| multiplication is done by columns of B and 
the result is stored in B» 

The following table summarizes the results of the program for 
various values of IFLAG. 

COMPUTE 

T 
A • B 

A • B 

T 
A • B 

A • B 

lAC is the number of rows of A 

JAB is the number of columns of 
(or B^) to be used, 
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STORE 


RESULT 


IN 




C 






C 






B 






B 




and C 


to be 


used. 


A and 


rows of B 



IBC is the number of columns of B (or bT) and 
C to be used* 

Called By ; LUNOBL, NUTPRE, OBLDRG, OBSERA, 
STACUL, STAPOS, XFORM 

DOMUD (TEST) 

Purpose ; 

This subroutine checks for errors. 

Description ; 

The Overflow and Divide Check Indicators are checked. If cither 

Is on, and the argument TEST Is not equal to ztra, TEST Is assumed to 

be a BCD word and the message "ERROR IN 'TEST'" Is printed. 

Called By ; DERIV, KEPLER, NUTPRE, OBLDRG, OBSERA, 
PRINT A, RECT, STAPOS, XFORM 

ECHREF 

Purpose ; 

This subroutine tests criteria for changing the reference body 
when the Encke Integrator Is used. 

Description : 

Criteria for changing the reference body based on regions of 
Influence are tested. 

Called By ; EITGRA 

Subroutines Called ; DDOT, EPHEM, KEPLER, SERVCE 

ECLIPS (J, K)^ 

Purpose : 

This subroutine checks for a change In the Illumination of the 
vehicle and prints out an appropriate message when the vehicle 
changes between the three states; umbra, penumbra, or sunlight* 

Description ; 

This first argument In the calling sequence Indicates whether 
the reference body Is the Earth ( J"»l) , the Moon (J-2), or some other 
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body (J>2). Thtt lecond argument indicates whether the vehicle la In 
the umbra (K-1) , penumbra (K-2) ^ or sunlight (K-3)« Values of the 
second argument are saved from one entry to the next and hence a 
simple comparison Is made for determining a change In vehicle status. 
If a change has occurred » a message Is printed containing the refer- 
ence body, the time of change and the entered state (umbra, penumbra, 
sunlight). 

Called By t DERIV 

EINT (lENT) 

Purpose ; 

This subroutine performs the Integration when the Encke method 
Is employed. 

Description ! 

The Gill modification of Rungc-Kutta Is used for short-term In- 
tegration and as a starting procedure for the Nordsleck long-term 
Integration. 

Called By ; DERIV 

Subroutines Called : ATIM, EITGRA, OBSERA 

EITGRA 

Purpose : 

This subroutine serves as the sub-main program governing calls 
to the Integration routines In the Encke method. 

Description : 

If not in powered flight, the program checks for a change of re- 
ference bodies. The delta of integration is determined by the distance 
of the vehicle from the reference body and integration is performed up 
to the next activity time by calling the integration routine. Fre- 
quent checks are also made on the magnitude of the deviations from the 
nominal orbit. If the magnitudes increase beyond specified limits, a 
rectification of the orbit is accomplished by calling RECT. 
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Called By : MAINA 

Subroutines Called ; ECHREF, EINT, ERSTRE, KEPLER, RECT 

EPHEM 

Purpose: 

This subroutine evaluates the position and velocity vectors of 
each of the seven bodies with respect to the reference body* 

Description : 

Tabular planetary positions are read from an ephemeris tape and 
Everett's interpolation formula is fitted to six points* It is eval- 
uated to find the position vector and its derivative is evaluated to 
find the velocity. 

Called By ; CCKREF, DERIV, ECHREF, STaCUL 

Subroutines Called ; SERVCE 

EKSTRE (IRC) 

Purpose : 

This subroutine saves or restores time, position, velocity, and 
acceleration of the vehicle at any designated time for the Encke in- 
tegrator* 

Description : 

Depending on the value of the argument (IRC) the information is 
saved (IRC <^ 1) or restored (IRC > 1). 

Called By ; ATIM, EITGRA, MAINA, OBSERA 

FIX (KTEMP, ITEMP, KNAME) 

Purpose : 

This subroutine unpacks a word into five separate words* 

Description : 

The argument KTEMP is assumed to contain a number of the form 
VVWWXXYYZZ* The word is unpacked and stored as follows: 

KNAME - W 

ITEMP(l) - ZZ 
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ITEMP(2) - YY 
ITEMP(3) - XX 

ITEMP(4) - VA^ 
Called By : OBSERA 

INPUTA 

Purpose ; 

This subroutine reeds all data necessary for one run. 

Degcription ? 

The subroutine Initializes necessary data and reads In sections 
desired. 

Called By ; EXECA 

Subroutines Called : XFORM 

KEPLER 

Purpose ; 

This subroutine computes the two^body position and velocity 
vectors. 

Description : 

A Newton Raphson scheme Is used to determine the differential ec- 
centric anomaly. After convergence » the tvorbody position end velo- 
city vectors are evalueted. The subroutine uses Herrlck*s method. 

Called By : ATIM, DERIV, ECHREF/ EITGRA, OBSERA, STACUL 

Subroutines Called : SERVCE, DOMUD 

LUNOBL(K) 

Purpose ; 

This subroutine computes the acceleration due to lunar oblateness* 
Optionally, It can compute the llbratlon and effect of the eerth*s 
J20 term. 

Description ; 

When K > 1, the llbratlon matrix is computed and then precessed 
and nutated. 

ISl 



Whan K - 2« th« Earth's J20 obl«t«n«s0 t«ra is QslcuUttdi 
When K « 3| both lunar and Karth oblatanaaa (J20 tarn) art 

computed. 

Called By t DERIV, OBD 

Subroutines Called ; DDOT, D>frML, NUTPRE, SERVCE 

MAINA 

Purpose t 

This subroutine handles the main flov of the prograa* 

Description : 

The routine calls subroutines which determine the next activity 
time, govern the ii:)tegration« compute the observations^ and print the 
chosen information. 

Called By t EXECA 

Subroutines Called ! CITGRA, CRSTRl, IITCRA, ERSTRE, OBSIRA, 

PFINIT, PRINTA, lECT, SERVCE, STACUL, 
TIMNGA 

MODELA(K) 

PurFOi^e i 

This subroutine computes the index of refraction for the tropo- 
sphere or ionosphere. 

Description : 

The index of refraction is computed with the tropospheriQ model 
when K "> 1 and the ionospheric model when R ■ 2» 

Called By I OBSERA 

MVP RAG 

Purpose ! 

This subroutine computes the ferturbations due to draf in the 
atmospheres of Mars or Venus* 

Description ! 

The coefficient of dreg is deurmlnad by interpolation txom |iv«n 
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tabltts AS A function of altitude and th« velocity of the vehicle* 
Celled By t DERIV 
Subroutines Called t DDOT, SERVCE 

NUTPRE(K) 

Purpose t 

This subroutine computes the expressions for determining the 
Greenwich hour angle » as well as the rotational matrices between co- 
ordinate systems* 

Description : 

When K ■■ If the expressions used In the nutation matrix and the 
llbratlon matrix are computed* 

When K - 2, the rotation matrix through the Greenwich hour 
angle Is computed* 

When K ■■ Sf the precession-nutation matrix Is computed* 

Called By ; LUNOBL, OBLDRG, PRINTA, STAPOS, XFORM 

gubroutlo»s Called t DMTML, DOMUD 

OBD 

Purpose: 

This subroutine computes observations from on«*board Instrxraen-* 
tatlon. 

Descrlptlon t 

The subroutine uses present position of vehicle with respect to 
reference bodies , landmarks, or ground stations to determine obser-* 
vat Ions. 

Called By t OBSERA 

Subroutines Called t DDOT, LUNOBL, SERVCE, STAPOS 

OBLDRG 

Purpose s 

This subroutine computes the oblateness and air drag perturbations 
due to the Earth and Its atmosphere* 
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Description ; 

OblatenesB Is computed In a flexible algorithm which can be used 
for zonal, tesseral^ or sectorial harmonics of any order. The drag 
perturbation Is computed from the Harrls-Prleater tables (at high alti- 
tudes) and the U. S, Standard Atmosphere (at low altitudes). 

Called By ; DERIV 

Subroutines Called ; DDOT, DENSTY, D>frML, DOMUD, NUTPRE, SERVCE 

OBSERA 

Purpose: 

This subroutine computes the observables as seen from a given 
ground station. 

Description : 

Given present vehicle location and ground station location, both 
referenced to Inertlal space, the observed values are computed. When 
specified, refraction corrections are Included. When In mode A, the 
total time in sight at the station is also computed. 

Called By ; MAINA 

Subroutines Called : ATIM, CINT, CRSTRE, DDOT, DMTML, DOMUD, 

EINT, ERSTRE, FIX, KEPLER, MODELA, OBD, 
SERVCE, STAPOS 

PFINIT 

Purpose : 

This subroutine performs initialization procedures at the entry 
of a burn period. 

Description : 

The subroutine determines coefficients of six Chebyshev polyno- 
mials which are valid for the first burn period. 

Called By: MAINA 
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PFLCHT 

Purpose : 

This routine computes the effect of a thrust perturbation In the 
En eke method* 

Description : 

The thrust acceleration vector at the Initial burn tlme^ the vehi- 
cle mass and mass rate are converted to a set of coefficients for 6 
Chebyshev polynomials. These coefficients are determined In PFINIT, 
The subroutine^ PFLGHT^ uses these coefficients to describe the powered 
flight trajectory as a function of time. Effectively^ the subroutine 
replaces KEPLER In Encke's method* Integration of the equations of 
motion continues exactly as If powered flight was not Involved* 

Called By : DERIV 

Subroutines Called; SERVCE 

PRINTA 

Purpose : 

This subroutine prints out current trajectory Information and 
observations* 

Description : 

Depending on whether the time entered corresponds to a print ^tlme^ 
the program checks each of the ten elements In the IPSEC array* If the 
value Is > O9 the corresponding section Is printed* 

To determine whether It Is a print tlme^ the subroutine first 
checks to see whether the present time Is within the print portion 
(DTPI) of a total print period (TAU). If not^ no printing Is done* 
If so. It next checks the value of the print Interval within DTPI 
(PRATE)* , If It Is negative. It always prints* If It Is positive, 
and It Is the first time Into the present print period. It prints, 
otherwise no printing Is done* When NDE - 1 or 4, printing occurs 
at each entry to the routine* 
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Called By : MAINA 

Subroutlneg Callad : DDOT, DOMUD, NUTPRE, SERVCE 

RECT 

Purpose ; 

This subroutine computes the parameters pertinent to a rectifi- 
cation In £ncke*s method. 

Description ; 

The two-body position and velocity vectors are equated to the 
true position and veolclty vectors. In addition, these components are 
saved. Perturbations in position and velocity are set equal to zero, 
and elements used by the KEPLER subroutine are computed. 

Called By ; EITGRA, MAINA 

Subroutines Called ; DDOT, DOMUD 

SERVCE (A, B, C, I) 

Purpose: 

This subroutine computes the cross product of two vectors as well 
as the magnitude, magnitude squared and magnitude cubed of a vector. 

Description ; 

When I « If the cross product of A and B is computed and 
the result stored in C(l), C(2), and C(3), and then continues as 
when I ■ 2. 

When 1-2, the magnitude cubed, magnitude, and magnitude 

squared is stored in C(4), C(5), and C(6), respectively. 

Called By ; CCHREF, DERIV, ECHREF, EPHEM, KEPLER, LUNOBL, 
MAINA, MVDRAG, OBD, OBLDRG, OBSERA, PFLGHT, 
PRINTA, STACUL 

STACUL 

Purpose : 

This subroutine determines the time of occultatlon (YCOM(23))t 
One of two types of occultatlon may be considered, vehicle or star 
occultatlon* 
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Description ; 

A Newton-Raphson Iteration Is used on the two-body solution to 
determine the time of occultatlon« When considering star occultatlon^ 
the occulting body Is the reference body, except In the Earth-Moon 
system, In which either the Moon or the Earth can do the occulting. 
Vehicle occultatlon occurs only In Moon reference « Up to three dlf-* 
ferent ground stations may be considered for occultatlon« 

Called By : MAINA 

Subroutines Called ; DDOT, DMTML, EPHEM, KEPLER, SERVCE, STAPOS 

STAPOS 

Purpose ; 

This subroutine computes the station position and velocity vectors. 

Description : 

Given the geodetic latitude, longitude, and altitude of the 
station, the coordinates of the station In the Greenwich coordinate 
system are computed. A rotation through the Greenwich hour angle then 
brings the coordinates Into the Inertlal system. The velocity of the 
station Is computed from the rotational speed of the earth. 

Called By ; OBD, OBSERA, STACUL 

Subroutines Called : D>rrML, DOMUD, NUTPRE 

TIMNGA 

Purpose : 

This subroutine determines the next activity time (I.e., the next 
time at which an observation Is to be computed or a printout of the 
trajectory Is to occur). 

Description : 

Logic Is set up for establishing an array of tines of Interest 
from which the earliest tine Is selected. Flags are set when the 
time selected Is the final tine or when e tine Is repeated. 
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Called By ; MAINA 

XFORM 

Purpose : 

This routine accepts the input information concerning the vehi- 
cle's initial conditions, and transforms them into the proper units 
and coordinate system* 

Description : 

After the initial conditions and associated flags are read in 
(Section 2 of the group I inputs), the program converts them into 
the position and velocity vectors in the units of the Earth radii, 
and Earth radii/hour. Then, if desired, the vectors are transformed 
into the base date system. 

Called By : INPUTA 

Subroutines Called: DMTML, DOMUD, NUTPRE 
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